{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An old FiPy solution to 1D BL\n",
    "I'm not really sure if the result is correct.\n",
    "The results is indeed incorrect. See the updated code a few cells below that might work slightly better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAEICAYAAAC5yopxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZo0lEQVR4nO3deXRc5Z3m8e9TJVmyjc1i2TR4iUxjCA5eIGLLQhsYCFsH+iSZZnUgdHyAMEn6NGlg+gw45AxNn2aabg4EjydDSDNhmQCHOMYNCQFPIJAOcsJmwGCwwcIJyKYx3m1Jv/mjSkUhZKusuqLqXj+fc4R0l6r7U6F6/N73vnVfRQRmZgC5WhdgZvXDgWBmJQ4EMytxIJhZiQPBzEocCGZW4kCwj42kcyX9vNZ12I7J4xB2T5ICmBIRy4fo+VuBFUBjRHQNxTEseW4h2KBIyte6BkueAyEDJF0h6S1J6yUtk3SCpCMlPSXpPUl/kHSzpGHF/X9VfOizkjZI+ktJF0h6os/zhqQDiz/fLulWSYskbQSOk3SapN9Lel/SKklzyx7ee4z3isc4pu8xJH1G0tOS1hW/f6Zs22JJ35P06+Lv9XNJLUPw8lkZB0LKSToYuAw4IiJGAV8AVgLdwF8DLcAxwAnApQARcWzx4TMiYo+IuKfCw50D/HdgFPAEsBGYDewFnAZcIunM4r69x9ireIyn+tS9D/AgcBMwBvgn4EFJY/oc70JgHDAMuLzCOm2QHAjp1w00AVMlNUbEyoh4LSKWRMRvIqIrIlYC/xP4syqP9dOI+HVE9ETElohYHBHPF5efA+7ahWOcBrwaEXcUa7wLeBn487J9fhgRr0TEZuD/AjOrrN8G4EBIuWKn4LeBucA7ku6WtL+kgyQtlPRHSe8D11FoLVRjVfmCpKMkPSapU9I64OJdOMb+wBt91r0BjC9b/mPZz5uAPXaxXttFDoQMiIg7I+JzwCeAAP4BuJXCv7hTImI08F8B7eRpNgIjehck/Ul/h+qzfCewAJgYEXsC88qOMdDlq9XFestNAt4a4HE2hBwIKSfpYEnHS2oCtgCbKZxGjALeBzZI+iRwSZ+Hvg0cULb8LPApSTMlNVNocQxkFPBuRGyRdCSFc/5enUBPn2OUWwQcJOkcSQ2S/hKYCiys4Lg2RBwI6dcEXA+sodDEHkehNXA5hTfoeuB/AX07DucCPypehfjPEfEKcC3wCPAqhU7DgVwKXCtpPXA1hfN8ACJiE4UOyF8Xj3F0+QMjYi1wOvA3wFrgb4HTI2JN5b+6Jc0Dk8ysxC0EMytxIJhZiQPBzEocCGZW0lCrA7e0tERra2utDm+221qyZMmaiBjb37aaBUJrayvt7e21OrzZbktS3xGiJT5lMLMSB4KZlTgQzKykZn0I/dm+fTsdHR1s2bKl1qVkSnNzMxMmTKCxsbHWpVidq6tA6OjoYNSoUbS2tiLt7IN5VqmIYO3atXR0dDB58uRal2N1rq5OGbZs2cKYMWMcBgmSxJgxY9zqsorUVSAADoMh4NfUKlV3gWBmteNAqMLtt9/O6tWrE3u+lStXcuedd5aW29vb+eY3v5nY85sNxIFQhcEEQlfXjucs6RsIbW1t3HTTTYOuz2xXORD62LhxI6eddhozZszg0EMP5Z577uHaa6/liCOO4NBDD2XOnDlEBPfeey/t7e2ce+65zJw5k82bN9Pa2sqaNYUb/rS3tzNr1iwA5s6dy5w5czjppJOYPXs2K1eu5POf/zyHH344hx9+OE8++SQAV155JY8//jgzZ87kxhtvZPHixZx++ukAvPvuu5x55plMnz6do48+mueee6703F/72teYNWsWBxxwgAPEqlJXlx3LffdnS3lx9fuJPufU/UdzzZ9/aqf7PPTQQ+y///48+OCDAKxbt44TTzyRq6++GoDzzz+fhQsX8uUvf5mbb76ZG264gba2tgGPvWTJEp544gmGDx/Opk2b+MUvfkFzczOvvvoqZ599Nu3t7Vx//fXccMMNLFxYuK3g4sWLS4+/5pprOOyww3jggQd49NFHmT17Ns888wwAL7/8Mo899hjr16/n4IMP5pJLLvGYAxsUtxD6mDZtGo888ghXXHEFjz/+OHvuuSePPfYYRx11FNOmTePRRx9l6dKlu/y8X/ziFxk+fDhQGID19a9/nWnTpvGVr3yFF198ccDHP/HEE5x//vkAHH/88axdu5Z169YBcNppp9HU1ERLSwvjxo3j7bff3uX6zKCCFoKk2yjcDPOdiDi0n+0C/gU4lcK98y+IiN9VW9hA/5IPlYMOOoglS5awaNEirrrqKk466SRuueUW2tvbmThxInPnzt3hNf2GhgZ6enoAPrLPyJEjSz/feOON7Lvvvjz77LP09PTQ3Nw8YF393fuy93JiU1NTaV0+n99pP4XZzlTSQrgdOHkn208BphS/5lCYDyC1Vq9ezYgRIzjvvPO4/PLL+d3vCtnW0tLChg0buPfee0v7jho1ivXr15eWW1tbWbJkCQD33XffDo+xbt069ttvP3K5HHfccQfd3d39Pl+5Y489lh//+MdA4VSipaWF0aNHV/fLmvUxYAshIn5VnNp7R84A/jUK/4T9RtJekvaLiD9UW9y2rh62dnV/aN1Hh9joQysHGoKj4n9yEg05kc/pQwN3nn/+eb7zne+Qy+VobGzk1ltv5YEHHmDatGm0trZyxBFHlPa94IILuPjiixk+fDhPPfUU11xzDRdddBHXXXcdRx111A5ruPTSS/nSl77ET37yE4477rhS62H69Ok0NDQwY8YMLrjgAg477LDSY+bOncuFF17I9OnTGTFiBD/60Y8G+E3Ndl1Ft2EvBsLCHZwyLASuj4gnisu/BK6IiI/c/UTSHAqtCCZNmvTpN9748H0aXnrpJQ455JDS8toNW3nrvc278OvsOiGaGnOMaMyzR3MDo5sbyeWyN7Kv72truy9JSyKi357wJK4y9Pfu6TdlImI+MB+gra1twCQaPbyR5sb8jp84PnqonT1p+e7dEXR1B109PWzZ3sO6Ldt5d9M2GnI5xo1uYszIYR7ya7udJAKhA5hYtjyBwrx9VWvM52jMfzwXQiKCDVu76Fy/ldXvbWbdpu1MGjPiYzu+WT1I4q99ATBbBUcD66rpP6jVTFKSGNXcyOSWkUzcewSbt3fzWucGtvXpw0gjz85llarksuNdwCygRVIHcA3QCBAR8yhM2nkqsJzCZccLB1tMc3Mza9eurelHoCWx98hhNDXmWLFmIyvWbOJPx46kIaUthd77IVRyadOskqsMZw+wPYBvJFHMhAkT6OjooLOzM4mnq9r2rh7+sGErnatyjNmjaeAH1KneOyaZDaSuhi43NjbW3V19fvjrFXz3vhe57i+mcc5Rk2pdjtmQSmc7+GP01WNaOeaAMfz9v73E2g1ba12O2ZByIAwglxPXnvEpNm3r5sZHXql1OWZDyoFQgSn7juLsIydyz9OrhnyglFktORAqdMmsAwGYt/i1GldiNnQcCBUav9dwzpw5nnuXdPD+lu21LsdsSDgQdsHsY1rZvL2b+5d01LoUsyHhQNgF0ybsyYwJe3LXb1fVuhSzIeFA2EVf+vQElr29nmV/7P++BWZp5kDYRaccuh85wc+eTe7262b1woGwi8aOauIzf9rCouervv+LWd1xIAzCCYeM4/U1G3lj7cZal2KWKAfCIMw6eBwAi5fVx4ewzJLiQBiEyS0j+cSYESxe9k6tSzFLlANhkD53YAu/XfEu3T2++YhlhwNhkNpa92bjtm5ffrRMcSAMUtsn9gFgyZv/UeNKzJLjQBikCXsPZ+yoJn73hgPBssOBMEiS+PSkvfm9WwiWIQ6EKhy07x68+e6mj8wuZZZWDoQqHDB2D3oCVr27qdalmCXCgVCFA8YW5mR8rdMjFi0bHAhVmNxSCITXHQiWEQ6EKoxqbmTsqCZe79xQ61LMEuFAqNIBLSNZscYtBMsGB0KVJuw9wnditsxwIFSpuTHH1q6eWpdhlggHQpWaGvJscyBYRlQUCJJOlrRM0nJJV/azfU9JP5P0rKSlkgY9A3TaDGvIORAsMwYMBEl54BbgFGAqcLakqX12+wbwYkTMoDB1/P+QNCzhWuvSsIYc27p76PHHoC0DKmkhHAksj4jXI2IbcDdwRp99AhglScAewLtAV6KV1qmmhsJLuK3brQRLv0oCYTxQPhFBR3FduZuBQ4DVwPPAtyLiI+8QSXMktUtq7+zMxu3HHAiWJZUEgvpZ17d9/AXgGWB/YCZws6TRH3lQxPyIaIuItrFjx+5ysfVoWG8guB/BMqCSQOgAJpYtT6DQEih3IXB/FCwHVgCfTKbE+jYsX3gJfenRsqCSQHgamCJpcrGj8CxgQZ993gROAJC0L3Aw8HqShdarpka3ECw7GgbaISK6JF0GPAzkgdsiYqmki4vb5wHfA26X9DyFU4wrImLNENZdN4bl84ADwbJhwEAAiIhFwKI+6+aV/bwaOCnZ0tKhtw/BN0mxLPBIxSo1uVPRMsSBUCVfZbAscSBUqXTK4HEIlgEOhCr1njJs3e5AsPRzIFTJIxUtSxwIVfJlR8sSB0KVfNnRssSBUCVfdrQscSBUyZcdLUscCFVyIFiWOBCq1JATOfnTjpYNDoQqSSrdRs0s7RwICRiW941WLRscCAkY1pD3KYNlggMhAU0NOY9DsExwICSgyXMzWEY4EBLgyVosKxwICSicMjgQLP0cCAlwC8GywoGQAI9DsKxwICTA4xAsKxwICWhqyPuyo2WCAyEB7kOwrHAgJMCBYFnhQEhAYz7Htu6+89+apY8DIQH5HPSEA8HSz4GQgLxEd48DwdKvokCQdLKkZZKWS7pyB/vMkvSMpKWS/l+yZda3fC7nQLBMGHCyV0l54BbgRKADeFrSgoh4sWyfvYDvAydHxJuSxg1VwfUon8OBYJlQSQvhSGB5RLweEduAu4Ez+uxzDnB/RLwJEBHvJFtmfcvlRLf7ECwDKgmE8cCqsuWO4rpyBwF7S1osaYmk2f09kaQ5ktoltXd2dg6u4jqUl+hxC8EyoJJAUD/r+v71NwCfBk4DvgD8N0kHfeRBEfMjoi0i2saOHbvLxdarvFsIlhED9iFQaBFMLFueAKzuZ581EbER2CjpV8AM4JVEqqxz+ZyIgJ6eIJfrLz/N0qGSFsLTwBRJkyUNA84CFvTZ56fA5yU1SBoBHAW8lGyp9SuvQgi4lWBpN2ALISK6JF0GPAzkgdsiYqmki4vb50XES5IeAp4DeoAfRMQLQ1l4PeltFXT3BI35GhdjVoVKThmIiEXAoj7r5vVZ/kfgH5MrLT3yxUDwaEVLO49UTEDplMFXGizlHAgJ6D1l6PEHHi3lHAgJaMi5U9GywYGQgN4WQpebCJZyDoQE9PYhOA8s7RwICcgXX0WfMljaORASkCu1EBwIlm4OhATkc77saNngQEhA3lcZLCMcCAlwC8GywoGQAI9UtKxwICQg5xaCZYQDIQGlcQjuQ7CUcyAkwH0IlhUOhAT448+WFQ6EBPQGQpenc7OUcyAkIOdbqFlGOBASkPf9ECwjHAgJ8IebLCscCAnwh5ssKxwICWjIFV5GX3a0tHMgJKCYB3Q5ECzlHAgJ8DgEywoHQgL84SbLCgdCAnJuIVhGOBAS4BaCZYUDIQH+cJNlRUWBIOlkScskLZd05U72O0JSt6QvJ1di/XMgWFYMGAiS8sAtwCnAVOBsSVN3sN8/UJglerfieypaVlTSQjgSWB4Rr0fENuBu4Ix+9vsvwH3AOwnWlwoeqWhZUUkgjAdWlS13FNeVSBoP/AXwoSni+5I0R1K7pPbOzs5drbVu+ZTBsqKSQFA/6/r+5f8zcEVEdO/siSJifkS0RUTb2LFjK62x7pWuMjgPLOUaKtinA5hYtjwBWN1nnzbgbhXeGC3AqZK6IuKBRKqsc/l8bwvBn3+2dKskEJ4GpkiaDLwFnAWcU75DREzu/VnS7cDC3SUMoHwcQo0LMavSgIEQEV2SLqNw9SAP3BYRSyVdXNy+036D3UHvh5s8UtHSrpIWAhGxCFjUZ12/QRARF1RfVrp4pKJlhUcqJsBXGSwrHAgJkITkUwZLPwdCQhpy8g1SLPUcCAnJSR6paKnnQEhIPif3IVjqORASkpf84SZLPQdCQnI5nzJY+jkQEpLPuYVg6edASIj7ECwLHAgJycuBYOnnQEhIoYVQ6yrMquNASEgu55GKln4OhIT4lMGywIGQkJyvMlgGOBAS0pAT3b6HmqWcAyEhOY9UtAxwICQk75GKlgEOhIR4pKJlgQMhITlfZbAMcCAkJJ+TxyFY6jkQEpLPiS5fZbCUcyAkJC+3ECz9HAgJ8acdLQscCAkpjFSsdRVm1XEgJCQvTwdv6edASIhPGSwLHAgJcSBYFlQUCJJOlrRM0nJJV/az/VxJzxW/npQ0I/lS65tHKloWDBgIkvLALcApwFTgbElT++y2AviziJgOfA+Yn3Sh9c4TtVgWVNJCOBJYHhGvR8Q24G7gjPIdIuLJiPiP4uJvgAnJlln/3EKwLKgkEMYDq8qWO4rrduQi4N/62yBpjqR2Se2dnZ2VV5kCvmOSZUElgaB+1vX7ly/pOAqBcEV/2yNifkS0RUTb2LFjK68yBTxRi2VBQwX7dAATy5YnAKv77iRpOvAD4JSIWJtMeenh2Z8tCyppITwNTJE0WdIw4CxgQfkOkiYB9wPnR8QryZdZ/3L+tKNlwIAthIjoknQZ8DCQB26LiKWSLi5unwdcDYwBvi8JoCsi2oau7PrjFoJlQSWnDETEImBRn3Xzyn7+K+Cvki0tXYblc2zr8kwtlm4eqZiQ5sY8Wx0IlnIOhIQ0NeTo7gm6PJ+bpZgDISFNjYWXcotbCZZiDoSENDXkAdi6vbvGlZgNngMhIc3FFoL7ESzNHAgJKbUQHAiWYg6EhDQ1FPsQfMpgKeZASEiTTxksAxwICWl2p6JlgAMhIW4hWBY4EBLS26noPgRLMwdCQno7Fd1CsDRzICSkudGXHS39HAgJ+aCF4FMGSy8HQkI+6ENwC8HSy4GQkA+uMriFYOnlQEhI6ZTBLQRLMQdCQiQxrCHnTkVLNQdCgpoach6HYKnmQEiQb6NmaedASFBTQ86dipZqDoQENbkPwVLOgZCgpoa8P+1oqeZASFBzo1sIlm4OhAQVWggOBEsvB0KCmhrdqWjp5kBIUGEcglsIll4OhAQVxiG4hWDpVVEgSDpZ0jJJyyVd2c92SbqpuP05SYcnX2r982VHS7sBZ3+WlAduAU4EOoCnJS2IiBfLdjsFmFL8Ogq4tfh9t9LUkGf9li7eem8zI4fla12O7Yb2HN6IpEE/vpLp4I8ElkfE6wCS7gbOAMoD4QzgXyMigN9I2kvSfhHxh0FXlkKfPbCFO3/7Jp+9/tFal2K7qRV/f2pVj68kEMYDq8qWO/jov/797TMe+FAgSJoDzAGYNGnSrtZa904+9E94+NvH8sSrnUStizEbhEoCob/2R9+/90r2ISLmA/MB2traMvmeOXDcHhw4bo9al2E2KJV0KnYAE8uWJwCrB7GPmdW5SgLhaWCKpMmShgFnAQv67LMAmF282nA0sG536z8wy4IBTxkiokvSZcDDQB64LSKWSrq4uH0esAg4FVgObAIuHLqSzWyoVNKHQEQsovCmL183r+znAL6RbGlm9nHzSEUzK3EgmFmJA8HMShwIZlbiQDCzEgeCmZU4EMysxIFgZiUOBDMrcSCYWYkDwcxKHAhmVqLC55JqcGCpE3ijgl1bgDVDXE616r3Geq8PXGNSKqnxExExtr8NNQuESklqj4i2WtexM/VeY73XB64xKdXW6FMGMytxIJhZSRoCYX6tC6hAvddY7/WBa0xKVTXWfR+CmX180tBCMLOPiQPBzErqNhAGmmC2FiRNlPSYpJckLZX0reL6fST9QtKrxe9717jOvKTfS1pYj/UVa9pL0r2SXi6+nsfUU52S/rr4//gFSXdJaq51fZJuk/SOpBfK1u2wJklXFd8/yyR9oZJj1GUglE0wewowFThb0tTaVgVAF/A3EXEIcDTwjWJdVwK/jIgpwC+Ly7X0LeClsuV6qw/gX4CHIuKTwAwK9dZFnZLGA98E2iLiUArTD5xVB/XdDpzcZ12/NRX/Ls8CPlV8zPeL76udi4i6+wKOAR4uW74KuKrWdfVT508pzIq9DNivuG4/YFkNa5pQ/MM4HlhYXFc39RVrGA2soNipXba+Lurkg7lK96EwVcFC4KR6qA9oBV4Y6DXr+56hMK/KMQM9f122ENjx5LF1Q1IrcBjw78C+UZypqvh9XO0q45+BvwV6ytbVU30ABwCdwA+LpzY/kDSSOqkzIt4CbgDepDBh8bqI+Hm91NfHjmoa1HuoXgOhoslja0XSHsB9wLcj4v1a19NL0unAOxGxpNa1DKABOBy4NSIOAzZSH6cxABTPw88AJgP7AyMlnVfbqnbZoN5D9RoIdTt5rKRGCmHw44i4v7j6bUn7FbfvB7xTo/I+C3xR0krgbuB4Sf+njurr1QF0RMS/F5fvpRAQ9VLnfwJWRERnRGwH7gc+U0f1ldtRTYN6D9VrIFQywezHTpKA/w28FBH/VLZpAfDV4s9fpdC38LGLiKsiYkJEtFJ4zR6NiPPqpb5eEfFHYJWkg4urTgBepH7qfBM4WtKI4v/zEyh0etZLfeV2VNMC4CxJTZImA1OA3w74bLXotKmw8+RU4BXgNeDval1PsabPUWh2PQc8U/w6FRhDoSPv1eL3feqg1ll80KlYj/XNBNqLr+UDwN71VCfwXeBl4AXgDqCp1vUBd1Ho09hOoQVw0c5qAv6u+P5ZBpxSyTE8dNnMSur1lMHMasCBYGYlDgQzK3EgmFmJA8HMShwIZlbiQDCzkv8PCcNkFBFhFbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer at 0x7f6bfee784e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAEICAYAAAC5yopxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZo0lEQVR4nO3deXRc5Z3m8e9TJVmyjc1i2TR4iUxjCA5eIGLLQhsYCFsH+iSZZnUgdHyAMEn6NGlg+gw45AxNn2aabg4EjydDSDNhmQCHOMYNCQFPIJAOcsJmwGCwwcIJyKYx3m1Jv/mjSkUhZKusuqLqXj+fc4R0l6r7U6F6/N73vnVfRQRmZgC5WhdgZvXDgWBmJQ4EMytxIJhZiQPBzEocCGZW4kCwj42kcyX9vNZ12I7J4xB2T5ICmBIRy4fo+VuBFUBjRHQNxTEseW4h2KBIyte6BkueAyEDJF0h6S1J6yUtk3SCpCMlPSXpPUl/kHSzpGHF/X9VfOizkjZI+ktJF0h6os/zhqQDiz/fLulWSYskbQSOk3SapN9Lel/SKklzyx7ee4z3isc4pu8xJH1G0tOS1hW/f6Zs22JJ35P06+Lv9XNJLUPw8lkZB0LKSToYuAw4IiJGAV8AVgLdwF8DLcAxwAnApQARcWzx4TMiYo+IuKfCw50D/HdgFPAEsBGYDewFnAZcIunM4r69x9ireIyn+tS9D/AgcBMwBvgn4EFJY/oc70JgHDAMuLzCOm2QHAjp1w00AVMlNUbEyoh4LSKWRMRvIqIrIlYC/xP4syqP9dOI+HVE9ETElohYHBHPF5efA+7ahWOcBrwaEXcUa7wLeBn487J9fhgRr0TEZuD/AjOrrN8G4EBIuWKn4LeBucA7ku6WtL+kgyQtlPRHSe8D11FoLVRjVfmCpKMkPSapU9I64OJdOMb+wBt91r0BjC9b/mPZz5uAPXaxXttFDoQMiIg7I+JzwCeAAP4BuJXCv7hTImI08F8B7eRpNgIjehck/Ul/h+qzfCewAJgYEXsC88qOMdDlq9XFestNAt4a4HE2hBwIKSfpYEnHS2oCtgCbKZxGjALeBzZI+iRwSZ+Hvg0cULb8LPApSTMlNVNocQxkFPBuRGyRdCSFc/5enUBPn2OUWwQcJOkcSQ2S/hKYCiys4Lg2RBwI6dcEXA+sodDEHkehNXA5hTfoeuB/AX07DucCPypehfjPEfEKcC3wCPAqhU7DgVwKXCtpPXA1hfN8ACJiE4UOyF8Xj3F0+QMjYi1wOvA3wFrgb4HTI2JN5b+6Jc0Dk8ysxC0EMytxIJhZiQPBzEocCGZW0lCrA7e0tERra2utDm+221qyZMmaiBjb37aaBUJrayvt7e21OrzZbktS3xGiJT5lMLMSB4KZlTgQzKykZn0I/dm+fTsdHR1s2bKl1qVkSnNzMxMmTKCxsbHWpVidq6tA6OjoYNSoUbS2tiLt7IN5VqmIYO3atXR0dDB58uRal2N1rq5OGbZs2cKYMWMcBgmSxJgxY9zqsorUVSAADoMh4NfUKlV3gWBmteNAqMLtt9/O6tWrE3u+lStXcuedd5aW29vb+eY3v5nY85sNxIFQhcEEQlfXjucs6RsIbW1t3HTTTYOuz2xXORD62LhxI6eddhozZszg0EMP5Z577uHaa6/liCOO4NBDD2XOnDlEBPfeey/t7e2ce+65zJw5k82bN9Pa2sqaNYUb/rS3tzNr1iwA5s6dy5w5czjppJOYPXs2K1eu5POf/zyHH344hx9+OE8++SQAV155JY8//jgzZ87kxhtvZPHixZx++ukAvPvuu5x55plMnz6do48+mueee6703F/72teYNWsWBxxwgAPEqlJXlx3LffdnS3lx9fuJPufU/UdzzZ9/aqf7PPTQQ+y///48+OCDAKxbt44TTzyRq6++GoDzzz+fhQsX8uUvf5mbb76ZG264gba2tgGPvWTJEp544gmGDx/Opk2b+MUvfkFzczOvvvoqZ599Nu3t7Vx//fXccMMNLFxYuK3g4sWLS4+/5pprOOyww3jggQd49NFHmT17Ns888wwAL7/8Mo899hjr16/n4IMP5pJLLvGYAxsUtxD6mDZtGo888ghXXHEFjz/+OHvuuSePPfYYRx11FNOmTePRRx9l6dKlu/y8X/ziFxk+fDhQGID19a9/nWnTpvGVr3yFF198ccDHP/HEE5x//vkAHH/88axdu5Z169YBcNppp9HU1ERLSwvjxo3j7bff3uX6zKCCFoKk2yjcDPOdiDi0n+0C/gU4lcK98y+IiN9VW9hA/5IPlYMOOoglS5awaNEirrrqKk466SRuueUW2tvbmThxInPnzt3hNf2GhgZ6enoAPrLPyJEjSz/feOON7Lvvvjz77LP09PTQ3Nw8YF393fuy93JiU1NTaV0+n99pP4XZzlTSQrgdOHkn208BphS/5lCYDyC1Vq9ezYgRIzjvvPO4/PLL+d3vCtnW0tLChg0buPfee0v7jho1ivXr15eWW1tbWbJkCQD33XffDo+xbt069ttvP3K5HHfccQfd3d39Pl+5Y489lh//+MdA4VSipaWF0aNHV/fLmvUxYAshIn5VnNp7R84A/jUK/4T9RtJekvaLiD9UW9y2rh62dnV/aN1Hh9joQysHGoKj4n9yEg05kc/pQwN3nn/+eb7zne+Qy+VobGzk1ltv5YEHHmDatGm0trZyxBFHlPa94IILuPjiixk+fDhPPfUU11xzDRdddBHXXXcdRx111A5ruPTSS/nSl77ET37yE4477rhS62H69Ok0NDQwY8YMLrjgAg477LDSY+bOncuFF17I9OnTGTFiBD/60Y8G+E3Ndl1Ft2EvBsLCHZwyLASuj4gnisu/BK6IiI/c/UTSHAqtCCZNmvTpN9748H0aXnrpJQ455JDS8toNW3nrvc278OvsOiGaGnOMaMyzR3MDo5sbyeWyN7Kv72truy9JSyKi357wJK4y9Pfu6TdlImI+MB+gra1twCQaPbyR5sb8jp84PnqonT1p+e7dEXR1B109PWzZ3sO6Ldt5d9M2GnI5xo1uYszIYR7ya7udJAKhA5hYtjyBwrx9VWvM52jMfzwXQiKCDVu76Fy/ldXvbWbdpu1MGjPiYzu+WT1I4q99ATBbBUcD66rpP6jVTFKSGNXcyOSWkUzcewSbt3fzWucGtvXpw0gjz85llarksuNdwCygRVIHcA3QCBAR8yhM2nkqsJzCZccLB1tMc3Mza9eurelHoCWx98hhNDXmWLFmIyvWbOJPx46kIaUthd77IVRyadOskqsMZw+wPYBvJFHMhAkT6OjooLOzM4mnq9r2rh7+sGErnatyjNmjaeAH1KneOyaZDaSuhi43NjbW3V19fvjrFXz3vhe57i+mcc5Rk2pdjtmQSmc7+GP01WNaOeaAMfz9v73E2g1ba12O2ZByIAwglxPXnvEpNm3r5sZHXql1OWZDyoFQgSn7juLsIydyz9OrhnyglFktORAqdMmsAwGYt/i1GldiNnQcCBUav9dwzpw5nnuXdPD+lu21LsdsSDgQdsHsY1rZvL2b+5d01LoUsyHhQNgF0ybsyYwJe3LXb1fVuhSzIeFA2EVf+vQElr29nmV/7P++BWZp5kDYRaccuh85wc+eTe7262b1woGwi8aOauIzf9rCouervv+LWd1xIAzCCYeM4/U1G3lj7cZal2KWKAfCIMw6eBwAi5fVx4ewzJLiQBiEyS0j+cSYESxe9k6tSzFLlANhkD53YAu/XfEu3T2++YhlhwNhkNpa92bjtm5ffrRMcSAMUtsn9gFgyZv/UeNKzJLjQBikCXsPZ+yoJn73hgPBssOBMEiS+PSkvfm9WwiWIQ6EKhy07x68+e6mj8wuZZZWDoQqHDB2D3oCVr27qdalmCXCgVCFA8YW5mR8rdMjFi0bHAhVmNxSCITXHQiWEQ6EKoxqbmTsqCZe79xQ61LMEuFAqNIBLSNZscYtBMsGB0KVJuw9wnditsxwIFSpuTHH1q6eWpdhlggHQpWaGvJscyBYRlQUCJJOlrRM0nJJV/azfU9JP5P0rKSlkgY9A3TaDGvIORAsMwYMBEl54BbgFGAqcLakqX12+wbwYkTMoDB1/P+QNCzhWuvSsIYc27p76PHHoC0DKmkhHAksj4jXI2IbcDdwRp99AhglScAewLtAV6KV1qmmhsJLuK3brQRLv0oCYTxQPhFBR3FduZuBQ4DVwPPAtyLiI+8QSXMktUtq7+zMxu3HHAiWJZUEgvpZ17d9/AXgGWB/YCZws6TRH3lQxPyIaIuItrFjx+5ysfVoWG8guB/BMqCSQOgAJpYtT6DQEih3IXB/FCwHVgCfTKbE+jYsX3gJfenRsqCSQHgamCJpcrGj8CxgQZ993gROAJC0L3Aw8HqShdarpka3ECw7GgbaISK6JF0GPAzkgdsiYqmki4vb5wHfA26X9DyFU4wrImLNENZdN4bl84ADwbJhwEAAiIhFwKI+6+aV/bwaOCnZ0tKhtw/BN0mxLPBIxSo1uVPRMsSBUCVfZbAscSBUqXTK4HEIlgEOhCr1njJs3e5AsPRzIFTJIxUtSxwIVfJlR8sSB0KVfNnRssSBUCVfdrQscSBUyZcdLUscCFVyIFiWOBCq1JATOfnTjpYNDoQqSSrdRs0s7RwICRiW941WLRscCAkY1pD3KYNlggMhAU0NOY9DsExwICSgyXMzWEY4EBLgyVosKxwICSicMjgQLP0cCAlwC8GywoGQAI9DsKxwICTA4xAsKxwICWhqyPuyo2WCAyEB7kOwrHAgJMCBYFnhQEhAYz7Htu6+89+apY8DIQH5HPSEA8HSz4GQgLxEd48DwdKvokCQdLKkZZKWS7pyB/vMkvSMpKWS/l+yZda3fC7nQLBMGHCyV0l54BbgRKADeFrSgoh4sWyfvYDvAydHxJuSxg1VwfUon8OBYJlQSQvhSGB5RLweEduAu4Ez+uxzDnB/RLwJEBHvJFtmfcvlRLf7ECwDKgmE8cCqsuWO4rpyBwF7S1osaYmk2f09kaQ5ktoltXd2dg6u4jqUl+hxC8EyoJJAUD/r+v71NwCfBk4DvgD8N0kHfeRBEfMjoi0i2saOHbvLxdarvFsIlhED9iFQaBFMLFueAKzuZ581EbER2CjpV8AM4JVEqqxz+ZyIgJ6eIJfrLz/N0qGSFsLTwBRJkyUNA84CFvTZ56fA5yU1SBoBHAW8lGyp9SuvQgi4lWBpN2ALISK6JF0GPAzkgdsiYqmki4vb50XES5IeAp4DeoAfRMQLQ1l4PeltFXT3BI35GhdjVoVKThmIiEXAoj7r5vVZ/kfgH5MrLT3yxUDwaEVLO49UTEDplMFXGizlHAgJ6D1l6PEHHi3lHAgJaMi5U9GywYGQgN4WQpebCJZyDoQE9PYhOA8s7RwICcgXX0WfMljaORASkCu1EBwIlm4OhATkc77saNngQEhA3lcZLCMcCAlwC8GywoGQAI9UtKxwICQg5xaCZYQDIQGlcQjuQ7CUcyAkwH0IlhUOhAT448+WFQ6EBPQGQpenc7OUcyAkIOdbqFlGOBASkPf9ECwjHAgJ8IebLCscCAnwh5ssKxwICWjIFV5GX3a0tHMgJKCYB3Q5ECzlHAgJ8DgEywoHQgL84SbLCgdCAnJuIVhGOBAS4BaCZYUDIQH+cJNlRUWBIOlkScskLZd05U72O0JSt6QvJ1di/XMgWFYMGAiS8sAtwCnAVOBsSVN3sN8/UJglerfieypaVlTSQjgSWB4Rr0fENuBu4Ix+9vsvwH3AOwnWlwoeqWhZUUkgjAdWlS13FNeVSBoP/AXwoSni+5I0R1K7pPbOzs5drbVu+ZTBsqKSQFA/6/r+5f8zcEVEdO/siSJifkS0RUTb2LFjK62x7pWuMjgPLOUaKtinA5hYtjwBWN1nnzbgbhXeGC3AqZK6IuKBRKqsc/l8bwvBn3+2dKskEJ4GpkiaDLwFnAWcU75DREzu/VnS7cDC3SUMoHwcQo0LMavSgIEQEV2SLqNw9SAP3BYRSyVdXNy+036D3UHvh5s8UtHSrpIWAhGxCFjUZ12/QRARF1RfVrp4pKJlhUcqJsBXGSwrHAgJkITkUwZLPwdCQhpy8g1SLPUcCAnJSR6paKnnQEhIPif3IVjqORASkpf84SZLPQdCQnI5nzJY+jkQEpLPuYVg6edASIj7ECwLHAgJycuBYOnnQEhIoYVQ6yrMquNASEgu55GKln4OhIT4lMGywIGQkJyvMlgGOBAS0pAT3b6HmqWcAyEhOY9UtAxwICQk75GKlgEOhIR4pKJlgQMhITlfZbAMcCAkJJ+TxyFY6jkQEpLPiS5fZbCUcyAkJC+3ECz9HAgJ8acdLQscCAkpjFSsdRVm1XEgJCQvTwdv6edASIhPGSwLHAgJcSBYFlQUCJJOlrRM0nJJV/az/VxJzxW/npQ0I/lS65tHKloWDBgIkvLALcApwFTgbElT++y2AviziJgOfA+Yn3Sh9c4TtVgWVNJCOBJYHhGvR8Q24G7gjPIdIuLJiPiP4uJvgAnJlln/3EKwLKgkEMYDq8qWO4rrduQi4N/62yBpjqR2Se2dnZ2VV5kCvmOSZUElgaB+1vX7ly/pOAqBcEV/2yNifkS0RUTb2LFjK68yBTxRi2VBQwX7dAATy5YnAKv77iRpOvAD4JSIWJtMeenh2Z8tCyppITwNTJE0WdIw4CxgQfkOkiYB9wPnR8QryZdZ/3L+tKNlwIAthIjoknQZ8DCQB26LiKWSLi5unwdcDYwBvi8JoCsi2oau7PrjFoJlQSWnDETEImBRn3Xzyn7+K+Cvki0tXYblc2zr8kwtlm4eqZiQ5sY8Wx0IlnIOhIQ0NeTo7gm6PJ+bpZgDISFNjYWXcotbCZZiDoSENDXkAdi6vbvGlZgNngMhIc3FFoL7ESzNHAgJKbUQHAiWYg6EhDQ1FPsQfMpgKeZASEiTTxksAxwICWl2p6JlgAMhIW4hWBY4EBLS26noPgRLMwdCQno7Fd1CsDRzICSkudGXHS39HAgJ+aCF4FMGSy8HQkI+6ENwC8HSy4GQkA+uMriFYOnlQEhI6ZTBLQRLMQdCQiQxrCHnTkVLNQdCgpoach6HYKnmQEiQb6NmaedASFBTQ86dipZqDoQENbkPwVLOgZCgpoa8P+1oqeZASFBzo1sIlm4OhAQVWggOBEsvB0KCmhrdqWjp5kBIUGEcglsIll4OhAQVxiG4hWDpVVEgSDpZ0jJJyyVd2c92SbqpuP05SYcnX2r982VHS7sBZ3+WlAduAU4EOoCnJS2IiBfLdjsFmFL8Ogq4tfh9t9LUkGf9li7eem8zI4fla12O7Yb2HN6IpEE/vpLp4I8ElkfE6wCS7gbOAMoD4QzgXyMigN9I2kvSfhHxh0FXlkKfPbCFO3/7Jp+9/tFal2K7qRV/f2pVj68kEMYDq8qWO/jov/797TMe+FAgSJoDzAGYNGnSrtZa904+9E94+NvH8sSrnUStizEbhEoCob/2R9+/90r2ISLmA/MB2traMvmeOXDcHhw4bo9al2E2KJV0KnYAE8uWJwCrB7GPmdW5SgLhaWCKpMmShgFnAQv67LMAmF282nA0sG536z8wy4IBTxkiokvSZcDDQB64LSKWSrq4uH0esAg4FVgObAIuHLqSzWyoVNKHQEQsovCmL183r+znAL6RbGlm9nHzSEUzK3EgmFmJA8HMShwIZlbiQDCzEgeCmZU4EMysxIFgZiUOBDMrcSCYWYkDwcxKHAhmVqLC55JqcGCpE3ijgl1bgDVDXE616r3Geq8PXGNSKqnxExExtr8NNQuESklqj4i2WtexM/VeY73XB64xKdXW6FMGMytxIJhZSRoCYX6tC6hAvddY7/WBa0xKVTXWfR+CmX180tBCMLOPiQPBzErqNhAGmmC2FiRNlPSYpJckLZX0reL6fST9QtKrxe9717jOvKTfS1pYj/UVa9pL0r2SXi6+nsfUU52S/rr4//gFSXdJaq51fZJuk/SOpBfK1u2wJklXFd8/yyR9oZJj1GUglE0wewowFThb0tTaVgVAF/A3EXEIcDTwjWJdVwK/jIgpwC+Ly7X0LeClsuV6qw/gX4CHIuKTwAwK9dZFnZLGA98E2iLiUArTD5xVB/XdDpzcZ12/NRX/Ls8CPlV8zPeL76udi4i6+wKOAR4uW74KuKrWdfVT508pzIq9DNivuG4/YFkNa5pQ/MM4HlhYXFc39RVrGA2soNipXba+Lurkg7lK96EwVcFC4KR6qA9oBV4Y6DXr+56hMK/KMQM9f122ENjx5LF1Q1IrcBjw78C+UZypqvh9XO0q45+BvwV6ytbVU30ABwCdwA+LpzY/kDSSOqkzIt4CbgDepDBh8bqI+Hm91NfHjmoa1HuoXgOhoslja0XSHsB9wLcj4v1a19NL0unAOxGxpNa1DKABOBy4NSIOAzZSH6cxABTPw88AJgP7AyMlnVfbqnbZoN5D9RoIdTt5rKRGCmHw44i4v7j6bUn7FbfvB7xTo/I+C3xR0krgbuB4Sf+njurr1QF0RMS/F5fvpRAQ9VLnfwJWRERnRGwH7gc+U0f1ldtRTYN6D9VrIFQywezHTpKA/w28FBH/VLZpAfDV4s9fpdC38LGLiKsiYkJEtFJ4zR6NiPPqpb5eEfFHYJWkg4urTgBepH7qfBM4WtKI4v/zEyh0etZLfeV2VNMC4CxJTZImA1OA3w74bLXotKmw8+RU4BXgNeDval1PsabPUWh2PQc8U/w6FRhDoSPv1eL3feqg1ll80KlYj/XNBNqLr+UDwN71VCfwXeBl4AXgDqCp1vUBd1Ho09hOoQVw0c5qAv6u+P5ZBpxSyTE8dNnMSur1lMHMasCBYGYlDgQzK3EgmFmJA8HMShwIZlbiQDCzkv8PCcNkFBFhFbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from fipy import *\n",
    "\n",
    "u = 1.e-3\n",
    "L = 100.\n",
    "nx = 200\n",
    "dt = 200.\n",
    "dx = L/nx\n",
    "muo = 0.002\n",
    "muw = 0.001\n",
    "mesh = Grid1D(dx = L/nx, nx = nx)\n",
    "x = mesh.cellCenters\n",
    "sw = CellVariable(mesh=mesh, name=\"saturation\", hasOld=True, value = 0.)\n",
    "sw.setValue(1,where = x<=dx)\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u\n",
    "*(sw**2./muw)/(sw**2./muw+(1-sw)**2./muo) * [[1]]) == 0\n",
    "\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "steps = 100\n",
    "viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)\n",
    "for step in range(steps):\n",
    "    sw.updateOld()\n",
    "    swres = 1.0e6\n",
    "    while swres > 1e-5:\n",
    "        swres = eq.sweep(dt = dt, var = sw)\n",
    "        print(swres)\n",
    "        \n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Updates solution on October 8th, 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fipy import *\n",
    "\n",
    "# relperm parameters\n",
    "swc = 0.1\n",
    "sor = 0.1\n",
    "krw0 = 0.3\n",
    "kro0 = 1.0\n",
    "nw = 2.0\n",
    "no = 2.0\n",
    "\n",
    "# domain and boundaries\n",
    "k = 1e-12 # m^2\n",
    "u = 1.e-3\n",
    "L = 100.\n",
    "nx = 100\n",
    "dt = 800.\n",
    "dx = L/nx\n",
    "\n",
    "# fluid properties\n",
    "muo = 0.002\n",
    "muw = 0.001\n",
    "\n",
    "# define the fractional flow functions\n",
    "def krw(sw):\n",
    "    res = krw0*((sw-swc)/(1-swc-sor))**nw\n",
    "    return res\n",
    "\n",
    "def dkrw(sw):\n",
    "    res = krw0*nw/(1-swc-sor)*((sw-swc)/(1-swc-sor))**(nw-1)\n",
    "    return res\n",
    "\n",
    "\n",
    "def kro(sw):\n",
    "    res = kro0*((1-sw-sor)/(1-swc-sor))**no\n",
    "    return res\n",
    "\n",
    "def dkro(sw):\n",
    "    res = -kro0*no/(1-swc-sor)*((1-sw-sor)/(1-swc-sor))**(no-1)\n",
    "    return res\n",
    "\n",
    "def fw(sw):\n",
    "    res = krw(sw)/muw/(krw(sw)/muw+kro(sw)/muo)\n",
    "    return res\n",
    "\n",
    "def dfw(sw):\n",
    "    res = (dkrw(sw)/muw*kro(sw)/muo-krw(sw)/muw*dkro(sw)/muo)/(krw(sw)/muw+kro(sw)/muo)**2\n",
    "    return res\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "sw_plot = np.linspace(swc, 1-sor, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the relative permeability and fractional flow curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5f3/8dcne+9Fwggj7CEQUBFURAVxIIoVd9VKtbXDb4d+tbZ+21+tbb/229o6iqNqHRT3AATrRClK2ENG2CEJZJO9zvX74zrBCIGcQJL7nJPP8/E4j7Nuct6G+ObOdd/3dYkxBqWUUr4vwOkASimlOocWulJK+QktdKWU8hNa6Eop5Se00JVSyk8EOfXBSUlJJjMz06mPV0opn7R69epiY0xyW+85VuiZmZnk5OQ49fFKKeWTRGTv8d7TIRellPITWuhKKeUntNCVUspPaKErpZSf0EJXSik/0W6hi8gzInJIRDYd530RkUdEJFdENojIuM6PqZRSqj2e7KE/C8w4wfsXAVnu2zzg8VOPpZRSqqPaLXRjzKdA6Qk2mQU8b6yVQJyI9OqsgMeoPAhL7oGmhi77CKWU8kWdMYaeAexv9TzP/doxRGSeiOSISE5RUdHJfdr+lfDF4/D+/Sf355VSyk91RqFLG6+1uWqGMWa+MSbbGJOdnNzmlavtGz4LzvgefPEEbHz15L6GUkr5oc4o9DygT6vnvYH8Tvi6x3fBr6HvmfD2D+Dgli79KKWU8hWdUehvAze6z3Y5A6gwxhR0wtc9vsBguOpZCImChTdA3eEu/TillPIFnpy2+DLwH2CIiOSJyK0icruI3O7eZDGwC8gFngS+12VpW4tOs6Veuhve+h7o2qhKqR6u3dkWjTHXtPO+Ab7faYk6IvMsO/yy7D5Y8Qic9SNHYiillDfw/StFz/y+PVD67wdg93Kn0yillGN8v9BFYNajkDgIXr0ZDnft8VillPJWvl/oAKHRcPUL0FADC2+EpnqnEymlVLfzj0IHSB4Csx+HvFWw5G6n0yilVLfzn0IHO5Y++S5Y/Q9Y/ZzTaZRSqlv5V6EDnHc/DDwPFv8U8nTNUqVUz+F/hR4QCFc+DdG94F83QNUhpxMppVS38L9CB4hIgLkvQm0ZLLwJmhudTqSUUl3OPwsdIG0UXPZX2LcClv3C6TRKKdXl2r1S1KeNvgry18LKRyF9LIyZ63QipZTqMv67h97igl9D5hR4+4dwYI3TaZRSqsv4f6EHBsFVz0FUKiy4zq54pJRSfsj/Cx0gMhGueQnqyu10u3olqVLKD/WMQgd7kPTyx2D/F/YcdZ1uVynlZ/z7oOjRRsyGwk2w/H8hbTRMvM3pREop1Wl6zh56i6n3weAZ8N49Ot2uUsqv9LxCDwiAK+ZDwgB45SYo2+t0IqWU6hQ9r9ABwmJh7svQ3AQLroX6KqcTKaXUKeuZhQ6QNAjmPAOHtsAb3wWXy+lESil1SnpuoQNknQ8X/ha2vgsfP+h0GqWUOiU96yyXtpxxh91L//SPkDwURs1xOpFSSp2Unr2HDnZN0ov/BH0nwVvfh7zVTidSSqmTooUOEBQCV/8TolLsQVJdaFop5YO00FtEJsE1C6ChypZ6Q43TiZRSqkO00FtLHQFXPAn56+zwi04PoJTyIVroRxs6E87/FWx+HT7+ndNplFLKY3qWS1vO+jEU58Inv4fEQTD6W04nUkqpdukeeltE4JL/g36T7dDLvpVOJ1JKqXZpoR9Py5kvsX3sQdLS3U4nUkqpE9JCP5GIBLh2Ibia4aWrobbc6URKKXVcWujtSRpk99RLd8Ir34bmRqcTKaVUmzwqdBGZISLbRCRXRO5p4/1YEXlHRNaLyGYRubnzozqo/9lwyZ9h10ew5Od6OqNSyiu1e5aLiAQCjwIXAHnAKhF52xizpdVm3we2GGMuFZFkYJuIvGiMaeiS1E4YdwOU7IDP/2LnUp/0A6cTKaXUN3iyhz4RyDXG7HIX9AJg1lHbGCBaRASIAkqBpk5N6g2mPQDDZ8Gy+2HLW06nUUqpb/Ck0DOA/a2e57lfa+1vwDAgH9gI/MgYc8wE4yIyT0RyRCSnqKjoJCM7KCAAZv8demfD6/Ng/yqnEyml1BGeFLq08drRg8jTgXVAOnAa8DcRiTnmDxkz3xiTbYzJTk5O7nBYrxAcbud8iU6Dl+fq6YxKKa/hSaHnAX1aPe+N3RNv7WbgdWPlAruBoZ0T0QtFJsF1r4FphhevgppSpxMppZRHhb4KyBKR/iISAswF3j5qm33ANAARSQWGALs6M6jXSRoEc1+C8r3wr+uhqd7pREqpHq7dQjfGNAF3AkuBr4CFxpjNInK7iNzu3uw3wCQR2Qh8ANxtjCnuqtBeo98kuPxx2Pu5zs6olHKcR5NzGWMWA4uPeu2JVo/zgQs7N5qPGDUHyvfBB/8DcX1h2i+dTqSU6qF0tsXOMPkuO/Sy/GGISYcJ33E6kVKqB9JC7wwiMPNhqDwIi38G0b1g6MVOp1JK9TA6l0tnCQyCOU9D+lh49RbY/6XTiZRSPYwWemcKibSzM8ak29kZi3OdTqSU6kG00DtbZBJc/xpIALxwhR2GUUqpbqCF3hUSBsB1C6G6CF66CuornU6klOoBtNC7SsZ4uOo5KNwEC2+EJv+ZeFIp5Z200LvS4Avhskdg54fw5h3gOma+MqWU6jR62mJXG3s9VB2yFx5FJsOM39nTHJVSqpNpoXeHyXfZUv/icYhKgSn/5XQipZQf0kLvDiIw/UGoKf56T33cDU6nUkr5GS307hIQALMes1PtvvNDiEiEoTOdTqWU8iN6ULQ7BYXAt56HXqfBqzfD3hVOJ1JK+REt9O4WGgXXvQKxveGluVC40elESik/oYXuhMgkuOFNW+7/vAJKdjqdSCnlB7TQnRLXx5a6aYbnL4eKA04nUkr5OC10JyUPtvO+1JbBP2dDdYnTiZRSPkwL3WnpY+HaBXaBjBev1HlflFInTQvdG2ROtvO+FGyAl6+BxjqnEymlfJAWurcYMgNmPwF7PrMLZDQ3Op1IKeVjtNC9yehvwcw/wrZF7sm8mp1OpJTyIXqlqLeZeBs0VMG/H4DgcLj0EZ3MSynlES10bzT5LmiogU//AMERMOMhLXWlVLu00L3V1HuhoRpWPmrXKp32S6cTKaW8nBa6txKB6b+FxhpY/rDdUz/7p06nUkp5MS10byYCF//JlvqHv7F76mfc4XQqpZSX0kL3di3T7jbWwHv3QFAoZN/idCqllBfS0xZ9QWAQXPkMZE2Hd++CNf90OpFSygtpofuKlrnUB54Hb/8A1r3sdCKllJfRQvclwWEw9yXofza89T3Y8IrTiZRSXsSjQheRGSKyTURyReSe42xzroisE5HNIvJJ58ZURwSHwzULoO8keOO7sPkNpxMppbxEu4UuIoHAo8BFwHDgGhEZftQ2ccBjwGXGmBHAVV2QVbUIiYBr/wW9J8Br34Gv3nU6kVLKC3iyhz4RyDXG7DLGNAALgFlHbXMt8LoxZh+AMeZQ58ZUx2hZyi59LLzybdi2xOlESimHeVLoGcD+Vs/z3K+1NhiIF5GPRWS1iNzY1hcSkXkikiMiOUVFRSeXWH0tLMYukJE2Cv51A2xd7HQipZSDPCn0tiYRMUc9DwLGAxcD04H7RWTwMX/ImPnGmGxjTHZycnKHw6o2hMXCDW/YUl94I2xd5HQipZRDPCn0PKBPq+e9gfw2tnnPGFNtjCkGPgXGdE5E1a7wOLjxTeg12pb6V+84nUgp5QBPCn0VkCUi/UUkBJgLvH3UNm8BU0QkSEQigNOBrzo3qjqhlj31XqfZMfUtR/8VKaX8XbuFboxpAu4ElmJLeqExZrOI3C4it7u3+Qp4D9gAfAk8ZYzZ1HWxVZtaSj19HLx6M2x5y+lESqluJMYcPRzePbKzs01OTo4jn+336g7Di3MgLwfmPA0jZjudSCnVSURktTEmu6339EpRf9Ry9kvvCfDqrXpFqVI9hBa6vwqNtqXebxK8fptO6KVUD6CF7s9Co+DahTBwKrx9J6x6yulESqkupIXu70IiYO7LMPgiWPQT+M+jTidSSnURLfSeIDjMTr077DJYeq9d0k4p5Xe00HuKoBCY8w8YdRV88Gv48Lfg0BlOSqmuoUvQ9SSBQTD773YZu0//AA1VMP1Bu3apUsrnaaH3NAGBcOlfISQKVj4G9Yfh0kfs60opn6aF3hMFBMCMhyA0xu6p11fBFU/aYRmllM/SQu+pROC8++xFSMt+YYdfvvVPe1aMUson6UHRnm7SD+yQS+4H8MKVUFfhdCKl1EnSQlcw/iY750vel/DcpVCli48o5Yu00JU18kqY+xIUbYN/zICyvU4nUkp1kBa6+trg6XDDm1BdBM9Mh4NbnE6klOoALXT1Tf3OhJuX2IuO/jED9n3hdCKllIe00NWxUkfArcsgIgmenwXblzmdSCnlAS101bb4fnDLUkgeDC/PhfULnE6klGqHFro6vqhkuOldyDwL3vgufP6Izv+ilBfTQlcnFhYD171ql7F7/3547x5wNTudSinVBr1SVLUvKBSufAai02Hlo1BZALPn22l5lVJeQwtdeSYgAGY8CDHpsOw+e/HR3BchIsHpZEopNx1yUR0z6U6Y8wwcyIFnZkD5fqcTKaXctNBVx428Eq5/HSoL4anzoXCj04mUUmihq5PVfwrc8p6dR/2ZGbDjfacTKdXjaaGrk5c6HL7zb0gYAC9dDauedjqRUj2aFro6NTHpdqqAQefDov+CpfeBy+V0KqV6JC10depCo+xMjRNug//8DRbeAA01TqdSqsfRQledIzAIZv4Rpv8Oti6C5y6BqkNOp1KqR9FCV51HBM78Hlz9gp1698nzoHCT06mU6jG00FXnG3YJ3LwYXE12XvVtS5xOpFSPoIWuukbGOLjtQ0gcBC9fA5//RSf2UqqLeVToIjJDRLaJSK6I3HOC7SaISLOIzOm8iMpntZwBM3wWvP9LeOv70FTvdCql/Fa7hS4igcCjwEXAcOAaERl+nO1+Dyzt7JDKh4VEwJx/wDn3wLoX7YIZ1cVOp1LKL3myhz4RyDXG7DLGNAALgFltbPcD4DVAT21Q3xQQAFP/G658GvLXwpNTdboApbqAJ4WeAbSegSnP/doRIpIBzAaeONEXEpF5IpIjIjlFRUUdzap83ag59mBpcyM8fSFsfsPpREr5FU8KXdp47eijW38G7jbGnHDlA2PMfGNMtjEmOzk52dOMyp9kjId5H0PqSHjl2/DBr3XBDKU6iSfzoecBfVo97w3kH7VNNrBARACSgJki0mSMebNTUir/Ep0G334XFv8Mlj9sz1W/8kkIi3U6mVI+zZM99FVAloj0F5EQYC7wdusNjDH9jTGZxphM4FXge1rm6oSCQuHSv8DFD8POD+DJaVC8w+lUSvm0dgvdGNME3Ik9e+UrYKExZrOI3C4it3d1QOXHRGDCd+DGt6G2zF5ZunWx06mU8lliHLrYIzs72+Tk5Djy2coLle+3k3rlr4UpP4Wp99q51pVS3yAiq40x2W29p1eKKu8Q1wdufg/G3QjL/xdeuBKqS5xOpZRP0UJX3iM4DC77K1z6COxdAfPPgQNrnE6llM/QQlfeZ/xNdnk7sJN7rX7O2TxK+QgtdOWdMsbBvE+g31nwzg/hze/pohlKtUMLXXmvyES4/jU4++ew7iV4ahoUbXc6lVJeSwtdebeAQDjvPlvsVYdg/rmwYaHTqZTySlroyjcMmga3L4deY+D12+CdH0FjrdOplPIqWujKd8Skw03vwOS7YPWz8NQFULLT6VRKeQ0tdOVbAoPg/Afg2lfgcB78/WxYv8DpVEp5BS105ZsGXwi3f2aHYN74Lrw+D+ornU6llKO00JXviu1th2DO/W/Y+IrdW89f63QqpRyjha58W0AgnHsPfHuRXa/0qQtgxV/B5XI6mVLdTgtd+Yd+k+wQzODpsOwX8OIcqCx0OpVSxzh4uI5DlXVd8rW10JX/iEiAq1+wc6zv/RweOxO+esfpVErR2Oxi6eZCbn12FZMe+pCnl+/uks/xZMUipXxHyxzrmVPs+er/uh7G3gAzHoLQKKfTqR4m91AVC3P28/qaPIqrGkiJDuW7Zw/gW9l92v/DJ0ELXfmn5CFw67/h49/BZ/8Hez6DK+ZDn4lOJ1N+rrKukUUbCnhldR6r95YRFCCcNzSFqyf04ZzByQQFdt3AiBa68l9BIXD+ryDrAnj9u3bmxrN/Zm+BwU6nU37E5TKs3FXCK6vzWLKpgLpGFwOTI7l35lBmj+1NcnRot+TQQlf+r98kuOMzWHI3fPJ72P4eXP4EpA53OpnycftKanh1TR6vrc7jQHkt0WFBXDGuN1eN781pfeIQkW7No4WueoawWJj9BAy9GN75sV0849x7YNKP7NWnSnnocF0jizcU8NqaPFbtKUMEJg9K4uczhjB9RBphwc4tnag/yapnGXYp9D0TFv0XfPBruyj15Y9D8mCnkykv1uwyLN9RxGtrDrBscyH1TS4GJEfys+lDmD02g/S4cKcjAlroqieKTIKrnoNNr8Hin8Lfp8B598MZd+jC1OoIYwyb8w/z5toDvLU+n6LKemLDg/lWdh+uGJfhyJBKe7TQVc8kAqPmQOZkOwSz7D7Y8hbM+ps9Q0b1WPnltby1Lp831uax/WAVwYHCuUNSuGJsBucNSyE0yHv/0RdjjCMfnJ2dbXJychz5bKW+wRi7aMZ7d0NDNZzzczjrx3omTA9SUdPIkk0FvLUun5W7SzAGxvWNY/a43lwyqhfxkSFORzxCRFYbY7Lbek/30JUSgTFXw8DzYMnP4MP/B5vftHvr6WOdTqe6SF1jMx9uPcSbaw/w8bYiGppd9E+K5EfTspg9NoN+iZFOR+wwLXSlWkQlw1XPwsg5sOgn8OQ0mHSnnc0x2DsOeqlT09TsYsXOEt5en8/STYVU1jeRHB3K9Wf04/Kx6YzKiPW6cfGO0EJX6mjDLrFj6+/fD5//xY6tX/wnuwye8jkulyFnbxnvrM9n8cYCSqobiAoNYsbINC4/LYMzByYSGOC7Jd6aFrpSbQmPg8v+CqOugnfvgheusHvu0x+E6FSn06l2GGPYeKCCdzcU8O76fPIr6ggLDmDasFQuHZ3OuUOSHT1fvKtooSt1Iv3PhjtW2Plglj8MO9630wmMvxkCdLJSb9JymuGijQUs2lDAvtIaggOFs7OSufuioUwblkpUqH9Xnp7lopSninfYvfU9y6H3BLjkz5A20ulUPZoxhq2FlSzaUMCijQXsLq4mMEA4a1ASl4zuxfThacRG+NfZSnqWi1KdISnLLnm34V+w9F675N3E2+xB0/A4p9P1GC174os3FrBkUyG7i6sJEJg0MIl5Zw9g+og0ErzoNMPu5FGhi8gM4C9AIPCUMeaho96/Drjb/bQKuMMYs74zgyrlFURgzFzIutCe3vjF3+0Vp+f/D4y5RodhukjLmPjijYUs2VTA3pIaAgOEMwckctuUAVw4IpWkqO6Z0dCbtTvkIiKBwHbgAiAPWAVcY4zZ0mqbScBXxpgyEbkIeMAYc/qJvq4OuSi/kL/OTh+Qtwp6T4SZf4T005xO5ReaXYbVe8tYsqmApZsKya+oIyhAmDQoiZkj07iwh+6Jn+qQy0Qg1xizy/3FFgCzgCOFboxZ0Wr7lUDvk4+rlA9JPw1uWQbrX4L3fwXzz4Xsm2HqLyAy0el0PqehycXKXSW8t7mQZZsPUlxVT0hQAGdnJXHXBYO5YHgqcRE9r8Q95UmhZwD7Wz3PA060930rsKStN0RkHjAPoG/fvh5GVMrLBQTA2Oth6CXw0YOw6inY+BqcezdMuM0utKGOq7q+iU+2F7F0cyEfbj1EZV0TESGBTB2awowRaUwdmuL3Z6d0Fk++S22dcd/mOI2ITMUW+uS23jfGzAfmgx1y8TCjUr4hPA5m/gGyb7EHTZfeC6uehum/hcEz7Pi7AqCosp4Ptx5k6eaDfJZbTEOTi/iIYGaMSGP6iDQmZyX55XniXc2TQs8DWq9o2hvIP3ojERkNPAVcZIwp6Zx4SvmglKFww+v2nPWl98LLc6H/OfaipB56mqMxhp1FVSzbcpB/bznI2v3lGAMZceFcf3o/LhyRSna/+C5db7Mn8KTQVwFZItIfOADMBa5tvYGI9AVeB24wxmzv9JRK+aKsC2DAuZDzD/j4QTvv+phrYeq9EJvhdLou19jsYvXeMj746iDvbznInpIaAEZlxHLX+YOZNiyF4b1ifHruFG/TbqEbY5pE5E5gKfa0xWeMMZtF5Hb3+08AvwQSgcfcfzlNxzsKq1SPEhgMp8+zc68vfxi+nA+bXoXTb4fJd/nd+etl1Q18sr2ID7Ye4pNthzhc10RIYABnDkzkO1MGMG1YCr1idaKzrqJXiirVncr22vPXNy6E8HiY8lN7cVKQb55DbYxh28FKPtpaxEdbD5GztxSXgaSoEKYOSWHasFQmZyXpQc1OdKLTFrXQlXJCwXp7muOujyC2r12wevTVPrFgdU1DE5/nlvDRtkN8vPUQ+RV1AAzrFcP5w2yJj86IJcBPZjD0NlroSnmrnR/Cvx+wBZ+YBVP/G4bP9qorTlsOaH68rYhPthfxxa5SGppdRIYEMjkrialDUjh3SAppsWFOR+0RtNCV8mbGwFfvwEe/haKtkDoSpt4HQy5y7FTHyrpGVuws4ZPtRXyyrYgD5bUADEqJ4pzByZw3NIUJmQmEBHnPPzw9hRa6Ur7A1WznhfnoQSjbDRnj7RkxA6d1ebE3u+xcKcu3F7F8RzFr9pXR5DJEhgRy1qAkzh2SwtmDk+gdH9GlOVT7tNCV8iXNjbD+ZfjkD1CxH9LHwTl3w+DpnVrsB8pr+WxHEZ/uKObz3GLKaxoBGJkRw9lZyUzJSmZ8v3jdC/cyOn2uUr4kMBjG3Qij59o5Ypb/CV6+GtJGwzk/hyEXn9QYe0VNI//ZVcxnucV8nlvC7uJqAFJjQjl/WCpTspKYPCiJRJ210GfpHrpS3q65ETYstOexl+6ElOEw5Scw/PITnhVT29BMzt5SVuwsYUVuMRsPVOAyEBESyBkDEjlrkC3wwalRenGPD9EhF6X8QXMTbH4DPv0jFG+D+Ew480447ToIiaChycX6vHJW5JawYmcxa/eV09DsIihAGNMnjsmDkpiclcSY3nE6jOLDtNCV8icuF2xbjOuzPxNwYBW1wXEsCr+Mh8vOpqAxAhEYkR7DpIFJnDkwkQmZCXphjx/RMXSl/EB9UzMb8ir4cncpK3clk7PvJ4xq2sx3m99lTuPzXBa8kPwhV5Ew9QfEZAxxOq5ygBa6Ul6qpqGJNXvL+XJ3CV/sLmXt/nIamlwADEmN5uoJfTljwGmM7f9DqNpByOePkLnpZch9AYbMhDPugMzJOm1vD6JDLkp5iaLKelbvLWXVnjJy9pSyKf8wzS5DgMDIjFgmZiYwsX8CEzITiD/e0muHC+wCGznPQG0ppI6yxT7ySgjWKzn9gY6hK+VlXC7DruIqVu8tI2dPGTl7y46cRhgaFMCYPnFk94tnYv8ExveLJzosuGMf0Fhrz4xZ+TgUfQWRyTD2Brs8XpyuFubLtNCVclh1fRPr88pZs7eM1XvLWLOvnIpaeyFPfEQw4/slMCEznuzMBEZmxBAa1Emr9RgDuz6GL/4OO5ba54OnQ/atMGgaBOiqQL5GD4oq1Y1a9r7X7Ctn7b5y1u4rY/vBSlzufafBqVHMHJXGuL7xjO8XT/+kyK47D1wEBk61t/L9sPpZWPM8bH/P7qmPv9muhxqV0jWfr7qV7qErdYoOVdaxfn8F6/eXsz6vnHX7y6msawIgOiyI0/rEMbZvPGP7xjGuTzyxER0cPulsTQ2w9V07zr5nOQQEQdZ0GHcDDLrAJ6bw7cl0D12pTlJR08im/Ao25FWwIa+c9fvLj8wHHhggDEmN5pLR6Yzra0t8QFKk980LHhQCI6+wt6LtsPafdu6YbYsgKhXGXGP32pOynE6qOkj30JU6jsN1jWw+cJhNByrYcKCCjXnlR9bFBOibEMFpfeIY0yeOMb1jGZEeS3iIj45JNzfCjmWw9gXYvhRMM/SeCKO/BSNmQ2SS0wmVmx4UVaodpdUNbDpQwab8Clvi+RXsbVXeGXHhjMqIZVTvWEb3jmVURixxEcc5ddDXVRbC+gX2LJlDm+2QzMBpttyHzIQQnULXSVroSrm5XIZ9pTVsKTjMlvzDR+4LD9cd2aZPQjgj02MZmRHL8PQYRmXEktRTZyAs3GTXP934Khw+ACFRMPRiOzHYwPP03HYHaKGrHqmitpFthZVsLTzM1sJKthYcZlthJdUNzYAd8x6UHMXw9BiG9YpmZEYsI3rFOn/Q0hu5XLD3c9jwL7u6Ul05hETbVZVGXG734LXcu4UWuvJrdY3N5B6qYvvBSrYftPfbCiuPLJsGEBsezJC0aIalRTM8PYbhvWLJSo0iLNhHx7yd1NwIuz6BLW/as2Vqy+ye++AZMHQmDDofwmKdTum3tNCVX6htaGZnURU7i6rYcbCKHYcq2XGwij0l1UfO8Q4OFAYmRzEkLZqhaTEMTYtmaK9o0mLCdM7vrtDcCLs/dZf7IqgpsWPumZPtePvgGRDfz+mUfkULXfkMYwyl1Q3sLKpml7u8dxZVk3uoiv1lNbT8uAYGCP0SIxicEs3gtGiGpEYzJC2KfomRBAfqXN+OcDVDXg5sWwzbltg52wFSRkDW+XZYpu8ZENRDj0d0Ei105XVqGprYU1zDnpJqdhdXs6uomt3FVewqrj6ytiXYeU36J0UyKCWKrJRoe58aRWZipC7S4O1Kdtpi3/4e7FsJrkYIjrR774Om2YJPHKizQXaQFrpyRFV9E3tLqtlbUuO+VbPH/bygou4b26bGhJKZGMnAlCgGJkcxMDmSgclRpMeFE+htF+aojquvhN3LYecHkPsBlO22r8f0tgXffwpkTtHhGQ/olaKqSzQ2uyisqGN/aQ37SmvYX1bDvtJa9pfWsL+0hpLqhm9snxQVQr/ESM4ckEj/pEj6J0fSPymSzMRIInVFHf8WGh0LiMYAAAnVSURBVG0PmA6daZ+X7ISdH9qpB3Lfhw0L7OuxfW259z0T+pwOiYNOakHsnkr30NVxVdc3UVBRS355HfnltRworyWvrJYDZbXkldVQeLjuyMFIsOPaGXHh9EkIp29CBH0SIshMjKRfYgT9EiN1GTTVNmOgaKvdg9+zHPZ8ZudyBwiPh94ToM9Ee+VqxngIjXI2r8N0D10do6q+icKKWgor6imoqOXg4ToKKuoorKgjv8IWeMv0ri0CBHrFhpMRH84ZAxLpHW8f94m35d0rNowgPSCpOkoEUobZ2+nz7DnvJbmw/wvI+xL2f2mnJbAb2zlmep0G6afZ+16j7W8ASvfQ/YnLZSiraaC4qoHiqnqKKu3t4OE6DrnvW563XFzTWnxEMKkxYWTEhdMrLoz0uHDSY8PpFWsfp8WG6Rkkyhm1ZfYMmgOrIX8dFKyDygL3mwIJAyB1OKS0uiUM8MuZI3UP3Uc1uwwVtY2UVjdQXtNAabW9lVQ3UFLVQGl1/ZHHxVX2cbPr2H+gw4IDSI0JIyU6lGHpMZwzJJnUmDB6xYZ9414vslFeKzwesi6wtxaVB6FgvS33gvVw6Ct7Lryx664SGAJJg+04fOJASBj49eOIRL88u8ajQheRGcBfgEDgKWPMQ0e9L+73ZwI1wLeNMWs6OatPqm9qprKuyX1r5HCtva+sa+JwXSPlNY1U1H59K69t5LC7xA/XNXK8X6AiQgJJjAohITKUtNgwRmbEkBwdSlJU6Dfuk6NDiQ4N0otqlP+JToXoC2HwhV+/1lgLxdttuR/aYu8LN9rpCkyr30pDYyG+L8T2cd96Q5z7cUy6XbIv0PemgGi30EUkEHgUuADIA1aJyNvGmC2tNrsIyHLfTgced997JZfL0Ohy0dRsaGx20dDkot59a2hy0eB+ra6xmdrGZuoam6lvdFHX1ExtQzM1Dfb1moYmaurt85rGZqrrm6iut+Vd3WAfNzafeEgrQCAmPJjY8GDiwoOJCQ+mT3w4CZEhxEWEEB8R/I3HiVGhJEaG6N60Um0JDodeY+ytteZGKN9nz64pyYXSnXYFp7I99mBsQ+WxXys8wc4PH5Vi7yOTICwOwuPc9/Hux7EQHOG+hUNQmGNn5niyhz4RyDXG7AIQkQXALKB1oc8Cnjd2QH6liMSJSC9jTMGxX+7UfLztEL95dwsGwIDBXl1o7+0whcuYb9y33BpdhqZmF22MSnRYWHAAESFBhAcHEhFib1FhQSRGRhAVGkRUWBCRoUFEhQYRHea+hQa7H9v7mPBgokODvG8BBKX8TWCwHWpJHAhceOz7teVQkQcV++3YfNUh9+2gvd//hZ3WoKHKs88LCrflHhBk120NCAIJ+PrxuJtg0p2d+p8InhV6BrC/1fM8jt37bmubDOAbhS4i84B5AH37ntzK49FhwQxNiwEBsV/TfW+fBwQIgSIEBsg3HgcGCEGBQnBAgL0PDCAowN6HBAUQGnT0fSBhwfY+PCSQsOBAwoICCAsOJDw4UEtYKX8S7t7zTht54u2aG6Guwv4DUFf+9X1THTTUQGONHfZpuXc12aEeV8vN/byL1nD1pNDbaq6j93E92QZjzHxgPtizXDz47GOM72cX1lVKqW4XGGyHXrx0BSdPBnrygD6tnvcG8k9iG6WUUl3Ik0JfBWSJSH8RCQHmAm8ftc3bwI1inQFUdMX4uVJKqeNrd8jFGNMkIncCS7GnLT5jjNksIre7338CWIw9ZTEXe9rizV0XWSmlVFs8Og/dGLMYW9qtX3ui1WMDfL9zoymllOoIvY5bKaX8hBa6Ukr5CS10pZTyE1roSinlJxybPldEioC9J/nHk4DiTozTWbw1F3hvNs3VMZqrY/wxVz9jTHJbbzhW6KdCRHKONx+wk7w1F3hvNs3VMZqrY3paLh1yUUopP6GFrpRSfsJXC32+0wGOw1tzgfdm01wdo7k6pkfl8skxdKWUUsfy1T10pZRSR9FCV0opP+HVhS4iM0Rkm4jkisg9bbw/VET+IyL1IvJTL8p1nYhscN9WiMiYtr6OA7lmuTOtE5EcEZnsDblabTdBRJpFZI435BKRc0Wkwv39Wiciv/SGXK2yrRORzSLyiTfkEpGftfpebXL/XSZ4Qa5YEXlHRNa7v1/dMhusB7niReQN9/+TX4pIO8slecAY45U37FS9O4EBQAiwHhh+1DYpwATgt8BPvSjXJCDe/fgi4AsvyRXF18dNRgNbvSFXq+0+xM7qOccbcgHnAu92x89VB3PFYdf07et+nuINuY7a/lLgQ2/IBdwL/N79OBkoBUK8INcfgV+5Hw8FPjjVz/XmPfQji1MbYxqAlsWpjzDGHDLGrAIavSzXCmNMmfvpSuwKTt6Qq8q4f3qASNpYJtCJXG4/AF4DDnVDpo7k6m6e5LoWeN0Ysw/s/wdekqu1a4CXvSSXAaJFRLA7NaVAkxfkGg58AGCM2QpkikjqqXyoNxf68RaedlpHc90KLOnSRJZHuURktohsBRYBt3hDLhHJAGYDT9B9PP17PNP9q/oSERnhJbkGA/Ei8rGIrBaRG70kFwAiEgHMwP4D7Q25/gYMwy6LuRH4kTHG5QW51gNXAIjIRKAfp7jz582F7tHC0w7wOJeITMUW+t1dmsj9cW281tZC3W8YY4YClwO/6fJUnuX6M3C3Maa5G/K08CTXGuy8GWOAvwJvdnkqz3IFAeOBi4HpwP0iMtgLcrW4FPjcGFPahXlaeJJrOrAOSAdOA/4mIjFekOsh7D/M67C/oa7lFH9z8GjFIod468LTHuUSkdHAU8BFxpgSb8nVwhjzqYgMFJEkY0xXTl7kSa5sYIH9jZgkYKaINBljurJA281ljDnc6vFiEXnMS75feUCxMaYaqBaRT4ExwHaHc7WYS/cMt4BnuW4GHnIPN+aKyG7smPWXTuZy/3zdDOAeDtrtvp28rj5ocQoHFYKAXUB/vj6oMOI42z5A9x0UbTcX0Be7vuokb/p+AYP4+qDoOOBAy3Nv+Ht0b/8s3XNQ1JPvV1qr79dEYJ83fL+wwwcfuLeNADYBI53O5d4uFjtGHdnVf4cd+H49Djzgfpzq/rlP8oJccbgPzgK3Ac+f6ud67R668WBxahFJA3KAGMAlIj/GHkk+fNwv3A25gF8CicBj7r3OJtPFM755mOtK4EYRaQRqgauN+6fJ4VzdzsNcc4A7RKQJ+/2a6w3fL2PMVyLyHrABcAFPGWM2OZ3LvelsYJmxvz10OQ9z/QZ4VkQ2YodC7jZd+1uWp7mGAc+LSDP2rKVbT/Vz9dJ/pZTyE958UFQppVQHaKErpZSf0EJXSik/oYWulFJ+QgtdKaX8hBa6Ukr5CS10pZTyE/8fp83zLvH1rpgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU9b3+8fcnIQlbgMSEnQSQRXYIIaCeWq1aUY/FrRVQQVyoVjkuP1s9ntau52q17altXRAtKm6gFautVrRUxVpFEgirLCEIJCBJIBvZJpn5/v5ItCkGmcAkz8zkfl1XruSZeZK5HWZuv/N9NnPOISIikS/G6wAiIhIaKnQRkSihQhcRiRIqdBGRKKFCFxGJEp28euCUlBQ3ePBgrx5eRCQi5eTklDjnUlu6z7NCHzx4MNnZ2V49vIhIRDKz3Ue7T1MuIiJRQoUuIhIlVOgiIlFChS4iEiVU6CIiUeKYhW5mi82syMw2HeV+M7PfmVmemW0ws4zQxxQRkWMJZoT+JDD9S+4/Hxje9DUfeOTEY4mISGsdcz9059wqMxv8JavMAJa4xvPwfmhmvcysn3Nuf4gyioiEtUDAUVFbT2l1PaXVPsqr6ymr8eFrCODzO+obAtT7G798fkdmehJnjGjx2KATEooDiwYAe5stFzTd9oVCN7P5NI7iSUtLC8FDi4i0vboGP7sPVrP3UNNXaQ0FpdXsPVTD/vIaymrqac2lJW786slhW+jWwm0t/qc55xYBiwAyMzN1ZQ0RCTvlNfVs2VfBlv0VbN5XzpZ9FeQVHaYh8K/K6hwXw6CkrgxK7kpGei+Su8bTq2s8vbrGkdT0vWeXODrHxRIXG0N8bAydYo242BjiYg2zlmrzxIWi0AuAQc2WBwL7QvB3RUTaXMnhOj7YeZAP8g/ywc6D7Cqp+vy+3okJjO7fg7NH9WZEn0TSkrsyMKkrKd3j26yUT0QoCv1V4BYzWwpMBco1fy4i4aq23s97O0p4P6+ED3YeZNuBSgC6J3Ri6pBkLp88kDH9ezCmf09SExM8Tts6xyx0M3seOBNIMbMC4IdAHIBzbiHwOnABkAdUA/PaKqyIyPGo8fl5d3sRr238lL9/fIAqn5/OcTFMGZzMjEn9Oe3kFMb270Gn2Mg+NCeYvVxmHeN+B9wcskQiIiFQ1+Bn5cdFvLZxP29vLaLa5ye5WzzfmNif88f2Y+rQZBI6xXodM6Q8O32uiEhb2H2wiuc+2sOL2QUcqvKR0j2eSyYN4MJx/cgakhzxo/Avo0IXkYjX4A+wcmsRz3y4m/d2lBAbY5wzqjdXTk3n9GEpxMaE3wbMtqBCF5GIVVXXwDMf7uaJ9z/h04pa+vbozO3njOCKKYPo27Oz1/HanQpdRCJOZW09Sz7YzePv5VNaXc9pJ5/Ej2eM4exTekf1lMqxqNBFJGKU19Tz5Puf8Id/5FNR28BZI1NZcPZwMtKSvI4WFlToIhL2anx+Hnsvn8dW5VNZ18C5o/uw4GvDGD+wl9fRwooKXUTClnOO1zbu5+evb6WwrIbzxvTh1rNHMLp/D6+jhSUVuoiEpU2F5fzkz1v46JNDjO7Xg//71gSmDj3J61hhTYUuImGl5HAdv35zG0vX7CWpazw/v3Qc38oc1GF2PTwRKnQRCQvOOV5eV8iPXt1Mtc/PdacPYcHZw+nZJc7raBFDhS4iniuqqOWelzfyt4+LmJyexH2XjWdY7+5ex4o4KnQR8UzzUXldQ4DvXziKeacP0fTKcVKhi4gnDlTUcs/yjazc2jgq/+Xl4xmaqlH5iVChi0i7e2vLAe58cT219X6NykNIhS4i7abeH+BXK7bx6Kp8xg7owe9mTtKoPIRU6CLSLvaX13DLc+vI2V3KVdPS+P6Fo+kcF13nI/eaCl1E2ty724u5fVkutfV+fjtzIjMmDvA6UlRSoYtIm/EHHA/8bTsPvp3HiN6JPHRlhnZHbEMqdBFpE4frGrj1+XWs3FrENycP5CczxtIlXlMsbUmFLiIht6+shmufXMOOosP8ZMYY5pw62OtIHYIKXURCav3eMq5fkk2tz8/ia6bw1RGpXkfqMFToIhIyr2/cz+3LcklNTODZ66cyok+i15E6FBW6iJww5xwPv7OTX67YxuT0JB69ejIp3RO8jtXhqNBF5IT4A47v/2kjz3+0lxkT+3PfZeO1f7lHVOgictx8DQFuX5bLaxv3c/NZJ3Pn10dipkP4vaJCF5HjUu1r4MZn1rJqezH/c8EobjhjqNeROjwVuoi0WnlNPdc9uYa1e0q577JxXDElzetIggpdRFqpuLKOuYs/YkdRJQ/OzuCCcf28jiRNVOgiErTCshquenw1+8treHyu9jEPNyp0EQlKQWk1Mxd9SHlNPc9cN5XMwcleR5IjqNBF5Jj2l9cw+7HVlNfU89z10xg3sKfXkaQFMcGsZGbTzWybmeWZ2d0t3N/TzP5sZuvNbLOZzQt9VBHxwoGKWmYt+pDSKh9PXzdVZR7GjlnoZhYLPAScD4wGZpnZ6CNWuxnY4pybAJwJ/NrM4kOcVUTaWVFlLbMe+5DiyjqevDaLiYN6eR1JvkQwI/QsIM85l++c8wFLgRlHrOOARGs8oqA7cAhoCGlSEWlXJYfrmP3Yaj4tr+XJa7OYnJ7kdSQ5hmAKfQCwt9lyQdNtzT0IjAL2ARuBW51zgSP/kJnNN7NsM8suLi4+zsgi0tYOVfm46vHVFJRWs/iaKUzRBtCIEEyht3Qcrzti+TwgF+gPTAQeNLMeX/gl5xY55zKdc5mpqdrdSSQcVdTWc9Xjq9lVUsUf5k5h2tCTvI4kQQqm0AuAQc2WB9I4Em9uHrDcNcoDdgGnhCaiiLSX2no/1z+VzY6iSh69ejKnD0vxOpK0QjCFvgYYbmZDmjZ0zgRePWKdPcDZAGbWBxgJ5IcyqIi0rQZ/gP96fh0f7TrEr745gTNH9vY6krTSMfdDd841mNktwAogFljsnNtsZjc23b8Q+CnwpJltpHGK5i7nXEkb5haREHLO8T8vb+LNLQf44UWjmTHxyM1kEgmCOrDIOfc68PoRty1s9vM+4OuhjSYi7eWXK7axLHsvC742jHmnD/E6jhynoA4sEpHo9Yd/7OLhd3YyKyuNO84d4XUcOQEqdJEO7E/rCvnpX7YwfUxffnbxWF2cIsKp0EU6qH/sKOHOF9czbWgyD8ycSGyMyjzSqdBFOqDtByq56ZkcTk7tzqI5mboGaJRQoYt0MMWVdcx7Yg2d42NZPG8KPTrHeR1JQkSFLtKB1Pj8XL8km4NVdfxhbiYDenXxOpKEkM6HLtJBBAKO25flsqGgjEevmsz4gTpzYrTRCF2kg7jvja28sflTvn/haL4+pq/XcaQNqNBFOoBnV+/m0VX5zDk1nWtPH+x1HGkjKnSRKPfejmLufWUzZ41M5d7/HK19zaOYCl0kiu0qqeLmZ9cyvHd3fj87g06xestHM/3rikSpitp6bliSTWyM8dicTLonaB+IaKd/YZEo5A84bluayyclVTx93VQGJXf1OpK0A43QRaLQr97cxt+3FvHDb4zh1JN1xaGOQoUuEmVeyS3kkXd2MntqGldPS/c6jrQjFbpIFNlQUMb3/riBrCHJ/OiiMV7HkXamQheJEkUVtcxfkkNK9wQeuTKD+E56e3c02igqEgV8DQG+8+xaymvqeemm0zipe4LXkcQDKnSRKPC/r20he3cpD86exOj+PbyOIx7RZzKRCPdSTgFPfbCbG74yhP8c39/rOOIhFbpIBNtUWM49L29k2tBk7pp+itdxxGMqdJEIVVrl48ZnckjuFs+DOqxf0By6SETyBxy3LsulqKKOF248lRRtBBVU6CIR6TdvbWfV9mJ+fuk4Jg7ShSqkkT6jiUSYFZs/5cG385g5ZRCzstK8jiNhRIUuEkF2lVRx5wvrGT+wJz/6ho4ElX+nQheJEDU+Pzc9k0NsrPHwlRl0jov1OpKEGc2hi0SIe1/ZxLYDlTxxzRQGJul0uPJFGqGLRIBla/bwYk4BC84axpkje3sdR8KUCl0kzG0qLOcHr2zmK8NTuPWcEV7HkTAWVKGb2XQz22ZmeWZ291HWOdPMcs1ss5m9G9qYIh1TeU0933l2Lcld43ngionExugCz3J0x5xDN7NY4CHgXKAAWGNmrzrntjRbpxfwMDDdObfHzPSZUOQEOee488X17CurYdm3T9UZFOWYghmhZwF5zrl855wPWArMOGKd2cBy59weAOdcUWhjinQ8j67K560tB7jnglFMTk/yOo5EgGAKfQCwt9lyQdNtzY0AkszsHTPLMbM5Lf0hM5tvZtlmll1cXHx8iUU6gNX5B/nlim1cOK4f804f7HUciRDBFHpLk3buiOVOwGTgQuA84Adm9oWtN865Rc65TOdcZmpqaqvDinQEJYfrWPD8OtKSu/KLy8ZhpnlzCU4w+6EXAIOaLQ8E9rWwTolzrgqoMrNVwARge0hSinQQ/oDjtqW5lNfU8+S8LBI7x3kdSSJIMCP0NcBwMxtiZvHATODVI9Z5BfiKmXUys67AVODj0EYViX4P/j2Pf+SV8JMZY3TlIWm1Y47QnXMNZnYLsAKIBRY75zab2Y1N9y90zn1sZm8AG4AA8LhzblNbBheJNu/nlfDAyu1cOmkA38ocdOxfEDmCOXfkdHj7yMzMdNnZ2Z48tki4Kaqo5YLfvUdS13heueV0usbrrBzSMjPLcc5ltnSfXjUiHmvwB1jw/Dqq6vw8f0OGylyOm145Ih77zd+2s3rXIX5zxQSG90n0Oo5EMJ3LRcRDb28r4qG3dzIraxCXTBrodRyJcCp0EY/sK6vhjmW5nNI3kR9epItVyIlToYt4oL5p3tzXENDFKiRkNIcu4oFfrdhGzu5Sfj9rEkNTu3sdR6KERugi7Wzlxwd4dFU+V01L46IJ/b2OI1FEhS7SjgpKq7njhfWM6d+D71842us4EmVU6CLtxNcQ4Jbn1uEPOB6arXlzCT3NoYu0k/ve2Eru3jIevjKDwSndvI4jUUgjdJF2sGLzp/zhH7uYe2o6F4zr53UciVIqdJE2tudgNXe+uJ7xA3tyz4WjvI4jUUyFLtKGauv9fOe5HAx4aHYGCZ00by5tR3PoIm3oZ69tYVNhBY/NyWRQclev40iU0whdpI28klvIMx/uYf4ZQzl3dB+v40gHoEIXaQN5RYf57+UbyUxP4rvnjfQ6jnQQKnSREKvx+bn52bV0jovl97MnERert5m0D82hi4TYD17ZxPaiSp6al0W/nl28jiMdiIYOIiH0QvZe/phTwIKzhnHGiFSv40gHo0IXCZHN+8r5wZ82cdrJJ3HrOSO8jiMdkApdJATKa+q56Zm19Ooax+9mTSI2xryOJB2Q5tBFTpBzjjtfXM++shqWfXsaKd0TvI4kHZRG6CIn6NFV+by15QD3XDCKyenJXseRDkyFLnICPth5kPvf2MqF4/sx7/TBXseRDk6FLnKciipqWfD8OgandOO+y8Zjpnlz8Zbm0EWOQ70/wM3PraWqroHnbphK9wS9lcR7ehWKHIf739jKmk9K+e3MiYzok+h1HBFAUy4irfbn9ft47L1dzDk1nRkTB3gdR+RzKnSRVtj6aQXf++MGMtOTdJFnCTsqdJEglVfX8+2nc0js3ImHr8wgvpPePhJeNIcuEoRAwHHbsnXsK6th6fxp9O7R2etIIl8Q1BDDzKab2TYzyzOzu79kvSlm5jezy0MXUcR7D6zcwdvbirn3ojE6eEjC1jEL3cxigYeA84HRwCwz+8LkYdN69wErQh1SxEtvbTnA71bu4PLJA7lqaprXcUSOKpgRehaQ55zLd875gKXAjBbWWwC8BBSFMJ+Ip3YWH+aOZbmMG9CTn108VgcPSVgLptAHAHubLRc03fY5MxsAXAIs/LI/ZGbzzSzbzLKLi4tbm1WkXVXWNm4EjesUw8KrJ9M5LtbrSCJfKphCb2lI4o5YfgC4yznn/7I/5Jxb5JzLdM5lpqbq5P8SvvwBx21Lc9lVUsWDsyYxoJeuPCThL5i9XAqAQc2WBwL7jlgnE1ja9HE0BbjAzBqcc38KSUqRdvbrN7excmsRP5kxhtOGpXgdRyQowRT6GmC4mQ0BCoGZwOzmKzjnhnz2s5k9CfxFZS6R6pXcQh5+ZyezstK4elq613FEgnbMQnfONZjZLTTuvRILLHbObTazG5vu/9J5c5FIsn5vGd/74wayBifz42+M0UZQiShBHVjknHsdeP2I21oscufcNSceS6T9FVXUMv/pbFK6J/DIVToSVCKPjhQVAWrr/cx/OofK2gZeuuk0TtJl5CQCqdClw3POcc/yjeTuLWPhVRmM6tfD60gix0WfKaXDe+TdnSxfV8jt54xg+th+XscROW4qdOnQ/rJhH/e/sY1vTOjPgq8N8zqOyAlRoUuHlbO7lDteWE9mehL3Xz6emBjt0SKRTYUuHdLug1XcsCSb/j07s2hOpg7rl6igQpcOp6zax7wn1xBwjifmZZHcLd7rSCIhoUKXDsXXEODbT+dQcKiGRVdnMiSlm9eRREJGuy1Kh+Gc4+7lG1i96xAPXDGRrCG6UIVEF43QpcP4zd92sHxt4+6JF08acOxfEIkwKnTpEJ7+cDe/W7mDb04eyH+drd0TJTqp0CXqvb5xP/e+somzT+nNzy8dpxNuSdRSoUtU+2DnQW5bmktGWhIPzs6gU6xe8hK99OqWqLV5Xznzl2QzOKUrf5ibSZd47Wsu0U2FLlFpz8Fq5i5eQ2LnTjx1bRa9umpfc4l+KnSJOsWVdVy9eDUNgQBLrsuiX09dD1Q6BhW6RJXy6nrmLv6IAxW1LL5mCsN6J3odSaTdqNAlalTW1jPniY/IKzrMo1dnkpGW5HUkkXalQpeoUFXXwLwn1rC5sJyHr8zgqyNSvY4k0u5U6BLxauv9XP9UNmv3lPLbmZM4Z3QfryOJeELncpGIVtfg59tP5/DhroP837cmcOF4XXFIOi6N0CVi1fsD3PLcOt7dXszPLxnHJZMGeh1JxFMqdIlI9f4Aty3N5a0tB/jxN8YwMyvN60gintOUi0ScugY/C55bx5tbDvA/F4xi7mmDvY4kEhZU6BJRauv93PhMDu9sK+ZHF43mmtOHeB1JJGyo0CViVPsauP6pbD7IP8jPLx3HLE2ziPwbFbpEhMraeuY9sYa1e0r59TcncGmGNoCKHEmFLmGvvLrxCNDNheX8flaGdk0UOQoVuoS14so65i5uPJz/kasmc64OGhI5KhW6hK384sPMfeIjSip9PDY3U4fzixxDUPuhm9l0M9tmZnlmdncL919pZhuavv5pZhNCH1U6knV7Srl84QdU1fl5fv40lblIEI45QjezWOAh4FygAFhjZq8657Y0W20X8FXnXKmZnQ8sAqa2RWCJfis/PsDNz62ld2Jnnro2iyEp3byOJBIRghmhZwF5zrl855wPWArMaL6Cc+6fzrnSpsUPAe2CIMdl6Ud7uGFJNsN7J/LSTaepzEVaIZg59AHA3mbLBXz56Ps64K8t3WFm84H5AGlp2odY/sU5x29X7uCBv+3gjBGpPHJlBt0StIlHpDWCecdYC7e5Flc0O4vGQv+Plu53zi2icTqGzMzMFv+GdDy19X7+e/lGXl5XyGUZA/nFZeOIi9VphkRaK5hCLwAGNVseCOw7ciUzGw88DpzvnDsYmngS7Q5U1DL/6RzW7y3jjnNHsOBrwzBraQwhIscSTKGvAYab2RCgEJgJzG6+gpmlAcuBq51z20OeUqJS7t4y5i/J5nBdAwuvmsz0sX29jiQS0Y5Z6M65BjO7BVgBxAKLnXObzezGpvsXAvcCJwEPN42uGpxzmW0XWyLd8rUF3L18I70TE1h+3Wmc0reH15FEIp45581UdmZmpsvOzvbkscU7/oDjvje2smhVPtOGJvPwlZNJ7hbvdSyRiGFmOUcbMGs3Amk3RZW13L4sl/fzDnL1tHTuvWi0Nn6KhJAKXdrF+3kl3Lo0l8raeu67bBxXTNFuqyKhpkKXNuUPNO5f/vu/72BoSjeevX4qI/smeh1LJCqp0KXNHKio5dal6/gw/xCXZQzkpxePoWu8XnIibUXvLmkT72wr4v+9sJ5qn59ffXMCl0/W2SBE2poKXULqcF0D//vaxzz/0R5G9OnO0tkZDO+jKRaR9qBCl5D5584SvvfHDRSW1fDtM4Zy+7kj6BwX63UskQ5DhS4nrNrXwP1vbOPJf37CkJRu/PHGU5mcnux1LJEOR4UuJ2TNJ4f47ovr+eRgNdecNpi7pp9Cl3iNykW8oEKX43Koysd9f93Ksuy9DEruwtL505g29CSvY4l0aCp0aZVAwLF0zV7uX7GVw7UNzD9jKLeePVznLhcJA3oXStA2FpTz/Vc2sX5vGVlDkvnZxWMZoT1YRMKGCl2O6VCVj9+8tZ1nVu/mpG4J/OaKCVw8cYDOWy4SZlToclTVvgYW/2MXC9/Np9rXwJxp6dzx9ZH07BLndTQRaYEKXb6gwR9gWfZeHvjbDoor6zh3dB/umj6SYb01vSISzlTo8jnnHCs2f8r9K7aRX1zF5PQkHrkyg8zB2qdcJBKo0AV/wPHXTft56O2dfLy/gmG9u/PYnEzOGdVb8+QiEUSF3oHV+wO8vK6Qhe/sJL+kiqGp3fjVNydw8cT+dNKFJ0Qijgq9A6rx+Xkhey+LVuVTWFbD6H49ePjKDM4b05fYGI3IRSKVCr0D2XOwmmdW72bZmr2U19STmZ7Ezy4Zy5kjUjW1IhIFVOhRLhBwvLujmKc/2M3b24qIMeO8MX2Ye+pgsoYkq8hFoogKPUodqKjlT+sKef6jPXxysJqU7gks+NpwZmel0bdnZ6/jiUgbUKFHkWpfA29uPsBLawt4P6+EgIPM9CTu+PpIpo/pS3wnbegUiWYq9AhX7w+wOv8QL68r5I1N+6ny+RnQqws3nzWMSyYNYGhqd68jikg7UaFHoBqfn1U7ilmx+VNWflxEeU09iQmd+M/x/bk0YwBTBicTo71VRDocFXqEKK6sY9X2Yt7c8invbi+mtj5Aj86dOGdUH74+pg9njuyty72JdHAq9DBV4/OzetdB3s8r4b0dJWz9tBKAvj06863MQZw3pi9ZQ5KJ0wFAItJEhR4mKmvryd1bRs7uUlbnHyJndyk+f4D42BgyByfxvekj+cqwVMb076HpFBFpkQrdA4GAY9fBKjYUlJH9SSk5u0vZdqAS58AMTunbg7mnpfMfw1PJGpysa3SKSFBU6G2stt7P9gOVbN5XwZZ9FWzeV87WTyup9vkB6J7QiUlpvZg+ti+T05OYOKgXiZ11vnERaT0VeggEAo6iyjrySw6TX1zFzuJ/fS8sq8G5xvUSEzoxqn8PvpU5iNH9ezBuQE9G9EnU+VNEJCSCKnQzmw78FogFHnfO/eKI+63p/guAauAa59zaEGf1RCDgKKupp7iyjuLKOg5U1FJYVkNhaQ0FZdUUltawr6wWnz/w+e90jY9laGo3MtKSuHzyQEb2SWRM/54MTOqi+W8RaTPHLHQziwUeAs4FCoA1Zvaqc25Ls9XOB4Y3fU0FHmn67innHP6Aw+cPUOPzU+3zU+VroKrOT7WvgWqfn8raBsqqfZTX1FNWXd/4vaae0iofxZV1lByuoyHgvvC3UxMTGJjUhbEDejJ9bD8GJHVhyEndOLl3N/r26KxzpIhIuwtmhJ4F5Dnn8gHMbCkwA2he6DOAJc45B3xoZr3MrJ9zbn+oA7+zrYifvfYxgYAj4BwBBwHncK7xQg31/gA+fwBfQ+N398UublGMQY8ucfTqEkfPrvEkd4vnlL6JpCYmkJqYQO/Ezp//3K9nZ+3zLSJhJ5hCHwDsbbZcwBdH3y2tMwD4t0I3s/nAfIC0tLTWZgUgsXMcI/skYgaxMUaMGWYQY0aMQVxsDHGxMSR0iiG+U8zny13jY+kaH0u3hE7/9j0xIY6eXeNITOik6RARiWjBFHpLLXfkuDeYdXDOLQIWAWRmZgY5dv53k9OTmJyedDy/KiIS1YI5zLAAGNRseSCw7zjWERGRNhRMoa8BhpvZEDOLB2YCrx6xzqvAHGs0DShvi/lzERE5umNOuTjnGszsFmAFjbstLnbObTazG5vuXwi8TuMui3k07rY4r+0ii4hIS4LaD9059zqNpd38toXNfnbAzaGNJiIiraFT9YmIRAkVuohIlFChi4hECRW6iEiUMBfssfGhfmCzYmD3cf56ClASwjihEq65IHyzKVfrKFfrRGOudOdcakt3eFboJ8LMsp1zmV7nOFK45oLwzaZcraNcrdPRcmnKRUQkSqjQRUSiRKQW+iKvAxxFuOaC8M2mXK2jXK3ToXJF5By6iIh8UaSO0EVE5AgqdBGRKBHWhW5m081sm5nlmdndLdx/ipl9YGZ1ZnZnGOW60sw2NH3908wmhEmuGU2Zcs0s28z+IxxyNVtvipn5zezycMhlZmeaWXnT85VrZveGQ65m2XLNbLOZvRsOuczsu82eq01N/5bJYZCrp5n92czWNz1f7XI22CByJZnZy03vyY/MbOwJP6hzLiy/aDxV705gKBAPrAdGH7FOb2AK8L/AnWGU6zQgqenn84HVYZKrO//abjIe2BoOuZqt93caz+p5eTjkAs4E/tIer6tW5upF4zV905qWe4dDriPWvwj4ezjkAu4B7mv6ORU4BMSHQa5fAj9s+vkUYOWJPm44j9A/vzi1c84HfHZx6s8554qcc2uA+jDL9U/nXGnT4oc0XsEpHHIddk2vHqAbLVwm0ItcTRYALwFF7ZCpNbnaWzC5ZgPLnXN7oPF9ECa5mpsFPB8muRyQaGZG46DmENAQBrlGAysBnHNbgcFm1udEHjScC/1oF572WmtzXQf8tU0TNQoql5ldYmZbgdeAa8Mhl5kNAC4BFtJ+gv13PLXpo/pfzWxMmOQaASSZ2TtmlmNmc8IkFwBm1hWYTuP/oMMh14PAKBovi7kRuNU5FwiDXOuBSwHMLAtI5wQHf+Fc6EFdeNoDQecys7NoLPS72jRR08O1cFtLF+p+2Tl3CnAx8NM2TxVcrgeAu5xz/nbI85lgcq2l8bwZEwXYL44AAAH9SURBVIDfA39q81TB5eoETAYuBM4DfmBmI8Ig12cuAt53zh1qwzyfCSbXeUAu0B+YCDxoZj3CINcvaPwfcy6Nn1DXcYKfHIK6YpFHwvXC00HlMrPxwOPA+c65g+GS6zPOuVVmdrKZpTjn2vLkRcHkygSWNn4iJgW4wMwanHNtWaDHzOWcq2j28+tm9nCYPF8FQIlzrgqoMrNVwARgu8e5PjOT9plugeByzQN+0TTdmGdmu2ics/7Iy1xNr695AE3TQbuavo5fW2+0OIGNCp2AfGAI/9qoMOYo6/6I9tsoesxcQBqN11c9LZyeL2AY/9oomgEUfrYcDv+OTes/SftsFA3m+erb7PnKAvaEw/NF4/TByqZ1uwKbgLFe52paryeNc9Td2vrfsBXP1yPAj5p+7tP0uk8Jg1y9aNo4C9wALDnRxw3bEboL4uLUZtYXyAZ6AAEzu43GLckVR/3D7ZALuBc4CXi4adTZ4Nr4jG9B5roMmGNm9UANcIVrejV5nKvdBZnrcuAmM2ug8fmaGQ7Pl3PuYzN7A9gABIDHnXObvM7VtOolwJuu8dNDmwsy10+BJ81sI41TIXe5tv2UFWyuUcASM/PTuNfSdSf6uDr0X0QkSoTzRlEREWkFFbqISJRQoYuIRAkVuohIlFChi4hECRW6iEiUUKGLiESJ/w8hvp1TJ/0RGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "krw_plot = [krw(sw) for sw in sw_plot]\n",
    "kro_plot = [kro(sw) for sw in sw_plot]\n",
    "fw_plot = [fw(sw) for sw in sw_plot]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(sw_plot, krw_plot, sw_plot, kro_plot)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(2)\n",
    "plt.plot(sw_plot, fw_plot)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the grid\n",
    "mesh = Grid1D(dx = L/nx, nx = nx)\n",
    "x = mesh.cellCenters\n",
    "\n",
    "# create the cell variables and boundary conditions\n",
    "sw = CellVariable(mesh=mesh, name=\"saturation\", hasOld=True, value = swc)\n",
    "# sw.setValue(1,where = x<=dx)\n",
    "sw.constrain(1-sor,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations\n",
    "$$\\varphi \\frac{\\partial S_w}{\\partial t}+u \\frac{\\partial f_w}{\\partial x}=0$$ or\n",
    "$$\\varphi \\frac{\\partial S_w}{\\partial t}+\\nabla.\\left( u \\frac{\\partial f_w}{\\partial S_w} S_w\\right)+ \\nabla. \\left( u f_w-u\\frac{\\partial f_w}{\\partial S_w} S_{w0} \\right)=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAEICAYAAABWCOFPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAd4klEQVR4nO3deXhV9b3v8fc380ACZGQImFACiowaBq0oYh2x2l7trdZKHVquWo+n59ZW7Xlu5Xie02uf41NbT61eb4+1ta32VnvUInVEKjhUgyLKHBkDQgZkJkCS7/1j79DsmJAN7GStJJ/X8+Qha+211/omsD/8fr+11m+ZuyMi0iIp6AJEJFwUCiISQ6EgIjEUCiISQ6EgIjEUCiISQ6Eg3cbMrjGzl4KuQ47OdJ1C32RmDpS7e1UX7b8UWA+kuntjVxxDuoZaCnJczCw56BqkaygUegEzu8PMtpjZHjNbbWbnmdkUM3vLzHaa2Sdm9nMzS4tu/3r0rR+Y2V4z+6qZXWdmi9vs181sZPT7x8zsITObb2b7gHPNbJaZvW9mu81ss5nNbfX2lmPsjB7jjLbHMLMzzexdM9sV/fPMVq8tNLN/NbM3oj/XS2ZW0AW/PmlDodDDmdlo4FZgsrvnABcCG4Am4J+AAuAM4DzgFgB3Pzv69gnu3s/d/xDn4b4G/BuQAywG9gGzgQHALOBmM/tSdNuWYwyIHuOtNnXnAc8DDwD5wE+A580sv83xrgeKgDTg9jjrlBOgUOj5moB0YIyZpbr7Bnf/2N2XuPvb7t7o7huA/wOcc4LHetbd33D3ZndvcPeF7v5hdHkZ8MQxHGMWsNbdH4/W+ASwCvhiq21+5e5r3P0A8P+AiSdYv8RBodDDRQcKvwPMBWrM7EkzG2Jmo8xsnpltM7PdwI+ItBpOxObWC2Y21cxeM7NaM9sF3HQMxxgCbGyzbiMwtNXytlbf7wf6HWO9chwUCr2Au//e3c8CTgIc+DHwEJH/ecvdPRf4AWBH2c0+IKtlwcwGtXeoNsu/B54Dhrl7f+DhVsfo7LTW1mi9rQ0HtnTyPuliCoUezsxGm9lMM0sHGoADRLoUOcBuYK+ZnQzc3Oat24ERrZY/AE41s4lmlkGk5dGZHGCHuzeY2RQiYwAtaoHmNsdobT4wysy+ZmYpZvZVYAwwL47jShdSKPR86cC9QB2R5nYRkVbB7UQ+pHuA/wu0HUycC/w6enbiv7v7GuAe4BVgLZGBxM7cAtxjZnuAHxLp9wPg7vuJDEq+ET3GtNZvdPd64FLgu0A98H3gUnevi/9Hl66gi5dEJIZaCiISQ6EgIjEUCiISQ6EgIjFSgjpwQUGBl5aWBnV4kT5ryZIlde5e2NHrgYVCaWkplZWVQR1epM8ys7ZXksZQ90FEYigURCSGQkFEYgQ2ptCew4cPU11dTUNDQ9Cl9CoZGRmUlJSQmpoadCnSA4QqFKqrq8nJyaG0tBSzo93QJ/Fyd+rr66murqasrCzocqQHCFX3oaGhgfz8fAVCApkZ+fn5an1J3EIVCoACoQvodyrHInShICLBUiicgMcee4ytW7cmbH8bNmzg97///ZHlyspKbrvttoTtXyQeCoUTcDyh0NjY8XNR2oZCRUUFDzzwwHHXJ3I8FApt7Nu3j1mzZjFhwgTGjh3LH/7wB+655x4mT57M2LFjmTNnDu7OU089RWVlJddccw0TJ07kwIEDlJaWUlcXmTiosrKSGTNmADB37lzmzJnDBRdcwOzZs9mwYQPTp0/ntNNO47TTTuPNN98E4M4772TRokVMnDiR+++/n4ULF3LppZcCsGPHDr70pS8xfvx4pk2bxrJly47s+4YbbmDGjBmMGDFCISInLFSnJFv7lz8vZ8XW3Qnd55ghudz9xVOPus0LL7zAkCFDeP755wHYtWsX559/Pj/84Q8BuPbaa5k3bx5XXnklP//5z7nvvvuoqKjo9NhLlixh8eLFZGZmsn//fl5++WUyMjJYu3YtV199NZWVldx7773cd999zJsXmaZw4cKFR95/9913M2nSJJ555hkWLFjA7NmzWbp0KQCrVq3itddeY8+ePYwePZqbb75Z1yTIcVNLoY1x48bxyiuvcMcdd7Bo0SL69+/Pa6+9xtSpUxk3bhwLFixg+fLlx7zfyy67jMzMTCBykda3vvUtxo0bx1e+8hVWrFjR6fsXL17MtddeC8DMmTOpr69n165dAMyaNYv09HQKCgooKipi+/btx1yfSItOWwpm9iiRCTZr3H1sO68b8DPgEiJz81/n7u+daGGd/Y/eVUaNGsWSJUuYP38+d911FxdccAEPPvgglZWVDBs2jLlz53Z4zj8lJYXm5maAz2yTnZ195Pv777+f4uJiPvjgA5qbm8nIyOi0rvbm0mw51Zienn5kXXJy8lHHLUQ6E09L4THgoqO8fjFQHv2aQ+R5Az3W1q1bycrK4utf/zq33347770XybeCggL27t3LU089dWTbnJwc9uzZc2S5tLSUJUuWAPD00093eIxdu3YxePBgkpKSePzxx2lqamp3f62dffbZ/O53vwMi3YqCggJyc3NP7IcVaUenLQV3fz36WPGOXA78xiP/lb1tZgPMbLC7f5KgGrvVhx9+yPe+9z2SkpJITU3loYce4plnnmHcuHGUlpYyefLkI9ted9113HTTTWRmZvLWW29x9913c+ONN/KjH/2IqVOndniMW265hSuuuII//vGPnHvuuUdaEePHjyclJYUJEyZw3XXXMWnSpCPvmTt3Ltdffz3jx48nKyuLX//61133S5A+La4p3qOhMK+D7sM84F53XxxdfhW4w90/M4OKmc0h0ppg+PDhp2/cGDvXw8qVKznllFOO/aeQTul3Ky3MbIm7dzg6noiBxvauoW03adz9EXevcPeKwsIOZ4MSkQAlIhSqgWGtlkuIPCdQRHqgRITCc8Bsi5gG7DqR8QQ9sSrx9DuVYxHPKckngBlAgZlVA3cDqQDu/jCRB4VeAlQROSV5/fEWk5GRQX19vW6fTqCW+RTiOe0pAvGdfbi6k9cd+HYiiikpKaG6upra2tpE7E6iWmZeEolHqC5zTk1N1exAIgHTZc4iEkOhICIxFAoiEiPUobB62x5+8tJqnVIT6UahDoVl1Tt5YEEVH21J7LwKItKxUIfCF04pJjnJeHH5tqBLEekzQh0KA7PTmFKap1AQ6UahDgWAC04tZm3NXtbV7g26FJE+oQeEwiAAXlyuKcZEukPoQ2HogEzGDe3PSyvUhRDpDqEPBYALTy3m/U07WbO9/anKRCRxekQofHXycPKz07jtifdpONwUdDkivVqPCIXCnHTu+8oEVm3bw71/WRV0OSK9Wo8IBYBzTy7ihs+X8dibG/jOk+9Ts1uPVhfpCqG6dbozd11yMv3Sk3n4r+t4ZWUNN5xVxo2fL6N/lp6GJJIocc3m3BUqKiq8svIzEz7HZX3dPn78l1W8sHwbOekpXDPtJK47s5RB/TW7kEhnOpvNuUeGQosVW3fz4GtV/OWjT0hOMi4aO5jZZ5xExUkDNZ2bSAd6dSi02LxjP796YwN/XLKZPQ2NjC7O4aopw/jypKEMyEpLyDFEeos+EQot9h9q5NmlW3ninU0sq95FWnIS548p5orThzK9vJDU5B4zrirSZfpUKLS2fOsu/lhZzbNLt/Dp/sPkZ6fxxQlDuHziECYOG6DuhfRZfTYUWhxqbGbh6hr+6/0tvLqqhkONzQzPy+LS8YO5dPwQThmco4CQPqXPh0Jruw4c5sWPtvHnZVt58+N6mpqdsoJsZo0bzMXjBjFmcK4CQno9hUIH6vce5IXl23h+2Se8va6eZofheVlcNHYQF55azKRhA0lKUkBI76NQiEP93oO8tGI7Ly7fxhtVdRxucgr6pXP+mCLOH1PMmZ8rICM1OegyRRJCoXCMdjcc5rVVNby0fDsLV9ew71ATWWnJnDWygC+MKebc0UUU5qQHXabIcessFHrUZc7dITcjlcsnDuXyiUM52NjEWx/X88rK7by6soaXVkQmepkwbAAzRxcx8+QiTh2Sq26G9CpqKcTJ3VnxyW4WrKzh1VU1fFC9E/fIHZwzRhUyY3QRZ5UX0D9T92FIuCWk+2BmFwE/A5KBX7r7vW1e7w/8FhhOpPVxn7v/6mj77Gmh0Fbd3oMsXF3LwtU1vL6mlt0NjSQnGacNH8A5owo5Z5RaERJOJxwKZpYMrAHOB6qBd4Gr3X1Fq21+APR39zvMrBBYDQxy90Md7benh0JrjU3NLN28MxISa2qOPKciPzuN6eUFTC8vZPqoAopydMOWBC8RYwpTgCp3Xxfd4ZPA5cCKVts4kGORk/z9gB1A43FX3cOkJCdRUZpHRWket184mto9B1lcVctfV9eyaG0dzyzdCsApg3M5OxoSFaUDdUZDQimelsKVwEXu/s3o8rXAVHe/tdU2OcBzwMlADvBVd3++nX3NAeYADB8+/PSNGzcm6ucIreZmZ/nW3SyqquX1NbUs2fgph5uc9JQkppTlMb28gLNGFurKSuk2iWgptPcvtW2SXAgsBWYCnwNeNrNF7h7zvDd3fwR4BCLdhziO3eMlJRnjSvozrqQ/t8wYyb6DjfxtfT2L1taxaG0dP5q/ClhFQb80Pj+ygLNGRloSmhtCghJPKFQDw1otlwBb22xzPXCvR5odVWa2nkir4Z2EVNmLZKenMPPkYmaeXAzAJ7sOsHhtHYur6nijqo5no12NkUX9OGtkAZ8fWcC0EXnkZOishnSPeLoPKUQGGs8DthAZaPyauy9vtc1DwHZ3n2tmxcB7wAR3r+tov71poDFRmpudVdv2sLiqlsVV9byzvp6Gw82kJBkThw3grPJIS2LCsAG6DVyOW6JOSV4C/JTIKclH3f3fzOwmAHd/2MyGAI8Bg4l0N+51998ebZ8Khc4dbGxiycZPj7QkPtyyC3fol57CtBH5nD0q0tUozc/SeITETZc59yI79x/izY8j4xGLq2rZvOMAACUDM5leXsjZ5QWcOVIXUMnRKRR6sQ11+1hUVcfra2p56+N69h6MXEA1aVjkAqoZo3UBlXyWQqGPONzUzPubdvL6mlpeX1vLsupdABT0S2fG6ELOHV3E9FEF5GrAss9TKPRRdXsP8vqaWl5bHbk+YteBw6QkGVNH5PGFU4o5f0wxJQOzgi5TAqBQEBqbmnlv004WrKrhlZXbqarZC8DYoblcOGYQl4wfzOcK+wVcpXQXhYJ8xvq6fby0fBsvLt/Ge5t2AnDyoBwumziEyyYMUQuil1MoyFFt29XAXz76hD9/sPVIQEwty+PK00uYNX4wWWmacqO3UShI3Dbv2M8z72/h6feq2VC/n5z0FL582lBmn3ESI4tygi5PEkShIMfM3Xl3w6c88c4mnv/wEw41NnPu6EJuOudzTB2RH3R5coIUCnJC6vYe5Hdvb+LxtzdQt/cQU8sit4dPLs0LujQ5TgoFSYgDh5p44p1NPPzXj6nZc5CLxw7iB5ecwrA8DUr2NJ2Fgu6qkbhkpiVzw1ll/PV75/I/zx/FwtW1XHD/6zy6eD1NzX3iLvg+Q6EgxyQzLZnbzivn1e+ew9QRedwzbwXX/uffqN97MOjSJEEUCnJchgzI5FfXTebHV4yjcuOnfPE/FrOsemfQZUkCKBTkuJkZX508nD/dfCZmxtWPvM27G3YEXZacIIWCnLCxQ/vz9M1nUtw/g288+g5vr6sPuiQ5AQoFSYhB/TN4cs40hgzIZM5vKtlUvz/okuQ4KRQkYYpyMnj0G5MxM/7Hb5dw4FBT0CXJcVAoSEINz8/ip1dNZNW23cx9bnnnb5DQUShIwp07uog5Z4/gD5WbeW/Tp0GXI8dIoSBd4h9mllOYk849f15Bsy5u6lEUCtIl+qWn8P0LR7N0806e/WBL0OXIMVAoSJe54rQSxg3tz09eXqPWQg+iUJAuk5RkfHN6GZt3HOCNjzt8LpCEjEJButSFpw5iQFYqT76zOehSJE4KBelSGanJ/LdJJby0YptumuohFArS5a6eMozDTc7T71UHXYrEQaEgXa68OIfTTxrIU0sUCj1BXKFgZheZ2WozqzKzOzvYZoaZLTWz5Wb218SWKT3dhacWs2b7Xmp2NwRdinSi01Aws2TgQeBiYAxwtZmNabPNAOAXwGXufirwlS6oVXqwadEJX9/SHZShF09LYQpQ5e7r3P0Q8CRweZttvgb8yd03Abh7TWLLlJ5uzOBcctJTeHud5lsIu3hCYSjQ+nxSdXRda6OAgWa20MyWmNns9nZkZnPMrNLMKmtra4+vYumRUpKTmFKWx9/UUgi9eEKhveeYt708LQU4HZgFXAj8LzMb9Zk3uT/i7hXuXlFYWHjMxUrPNm1EPuvq9rFd4wqhFk8oVAPDWi2XAFvb2eYFd9/n7nXA68CExJQovUXLuIJmZgq3eELhXaDczMrMLA24CniuzTbPAtPNLMXMsoCpwMrElio93ZghueRkpCgUQq7Tp4e6e6OZ3Qq8CCQDj7r7cjO7Kfr6w+6+0sxeAJYBzcAv3f2jrixcep7kJGNqWZ4GG0MurkcKu/t8YH6bdQ+3Wf534N8TV5r0RmOH9ufVVTU0HG4iIzU56HKkHbqiUbpVWUE27rBRE7uGlkJButWIgn4ArK/bF3Al0hGFgnSr0oLIA2kVCuGlUJBulZORSkG/dNbX7Q26FOmAQkG63YiCbLUUQkyhIN2urCCb9XUaaAwrhYJ0u9KCbOr2HmR3w+GgS5F2KBSk25UVZAOwQV2IUFIoSLcbURgJBY0rhJNCQbrd8LwszBQKYaVQkG6XkZrMkP6ZCoWQUihIIEYU6rRkWCkUJBCl+QqFsFIoSCAG9c9gT0MjDYebgi5F2lAoSCDystMA2LHvUMCVSFsKBQmEQiG8FAoSiJZQqFcohI5CQQLREgqfKhRCR6EggchXSyG0FAoSiNyMVJKTjB379Hj6sFEoSCCSkoyBWans2Kc7JcNGoSCByctOU0shhBQKEphIKGhMIWwUChIYhUI4KRQkMAqFcFIoSGDystPZeeAwTc1tH2IuQVIoSGDyslJxh5371VoIk7hCwcwuMrPVZlZlZnceZbvJZtZkZlcmrkTprfL6pQO6/yFsOg0FM0sGHgQuBsYAV5vZmA62+zGRp1OLdEpXNYZTPC2FKUCVu69z90PAk8Dl7Wz3D8DTQE0C65NebGCW7n8Io3hCYSiwudVydXTdEWY2FPgyEPN4+rbMbI6ZVZpZZW1t7bHWKr1Mfj+1FMIonlCwdta1HS7+KXCHux91Gh13f8TdK9y9orCwMN4apZdqaSloTCFcUuLYphoY1mq5BNjaZpsK4EkzAygALjGzRnd/JiFVSq+UlpJETkaKQiFk4gmFd4FyMysDtgBXAV9rvYG7l7V8b2aPAfMUCBIPXcAUPp2Ggrs3mtmtRM4qJAOPuvtyM7sp+vpRxxFEjkahED7xtBRw9/nA/Dbr2g0Dd7/uxMuSviI/O40tOxuCLkNa0RWNEqiBWWk6JRkyCgUJVF6/SPfBXfc/hIVCQQI1MCuNQ03N7D+kh8KEhUJBApWbkQrAnobGgCuRFgoFCVRuZmSse3eD5moMC4WCBCrnSEtBoRAWCgUJVG5GtKVwQN2HsFAoSKBaWgrqPoSHQkEC9fcxBbUUwkKhIIFqOfuw+4BaCmGhUJBAZaQmk5acpO5DiCgUJHC5mSm6TiFEFAoSuNyMVHUfQkShIIHLyVBLIUwUChK43MxUjSmEiEJBAqfuQ7goFCRw6j6Ei0JBAqfuQ7goFCRwuRkpNBxu5lBjc9ClCAoFCQHdKRkuCgUJnO5/CBeFggQuJ133P4SJQkECl5upKdnCRKEggdOUbOGiUJDA5ej26VBRKEjgWqZkU/chHBQKErjstBSSTN2HsIgrFMzsIjNbbWZVZnZnO69fY2bLol9vmtmExJcqvVVSktEvPUXdh5DoNBTMLBl4ELgYGANcbWZj2my2HjjH3ccD/wo8kuhCpXfLzUxV9yEk4mkpTAGq3H2dux8CngQub72Bu7/p7p9GF98GShJbpvR2uRm6/yEs4gmFocDmVsvV0XUduRH4S3svmNkcM6s0s8ra2tr4q5ReLycjRc9+CIl4QsHaWdfuI4LN7FwioXBHe6+7+yPuXuHuFYWFhfFXKb2e7pQMj5Q4tqkGhrVaLgG2tt3IzMYDvwQudvf6xJQnfUVuhsYUwiKelsK7QLmZlZlZGnAV8FzrDcxsOPAn4Fp3X5P4MqW3i3Qf1FIIg05bCu7eaGa3Ai8CycCj7r7czG6Kvv4w8EMgH/iFmQE0untF15UtvU1uZip7DzXS3OwkJbXXY5XuEk/3AXefD8xvs+7hVt9/E/hmYkuTvqR/ZirusOvAYQZmpwVdTp+mKxolFErzswBYV7cv4EpEoSChUF6UA8Da7XsCrkQUChIKJQMzyUhNYm3N3qBL6fMUChIKSUnGyKJ+rFFLIXAKBQmNUUU5rN2ulkLQFAoSGuXFOWzb3aArGwOmUJDQKC/qB6DWQsAUChIao4ojZyCqajSuECSFgoRGyxmINWopBEqhIKHRcgZCpyWDpVCQUCkvytEFTAGL694Hke5SXtyP/3p/Cz95eQ3TywvITE0OuqQe5dQhuURvSjxuCgUJlYvHDmbByhr+Y8FaHnh1bdDl9Djr//clJ7wPhYKESllBNk/dfCY79h1i6eZPaWxqd5Iv6UIKBQmlvOw0Zp5cHHQZfZIGGkUkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkRlyhYGYXmdlqM6syszvbed3M7IHo68vM7LTElyoi3aHTUDCzZOBB4GJgDHC1mY1ps9nFQHn0aw7wUILrFJFuEk9LYQpQ5e7r3P0Q8CRweZttLgd+4xFvAwPMbHCCaxWRbhBPKAwFNrdaro6uO9ZtMLM5ZlZpZpW1tbXHWquIdIN4QqG9WSDbzpEVzza4+yPuXuHuFYWFhfHUJyLdLJ5QqAaGtVouAbYexzYi0gPEEwrvAuVmVmZmacBVwHNttnkOmB09CzEN2OXunyS4VhHpBp1O3OrujWZ2K/AikAw86u7Lzeym6OsPA/OBS4AqYD9wfdeVLCJdKa7ZnN19PpEPfut1D7f63oFvJ7Y0EQmCrmgUkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJYZF7mQI4sFktsLGTzQqAum4o50SoxsRQjYkRT40nuXuHsxwFFgrxMLNKd68Iuo6jUY2JoRoTIxE1qvsgIjEUCiISI+yh8EjQBcRBNSaGakyME64x1GMKItL9wt5SEJFuplAQkRihDYXOHmobBDMbZmavmdlKM1tuZv8YXZ9nZi+b2dronwMDrjPZzN43s3lhrC9a0wAze8rMVkV/n2eErU4z+6fo3/NHZvaEmWUEXaOZPWpmNWb2Uat1HdZkZndFP0OrzezCeI4RylCI86G2QWgEvuvupwDTgG9H67oTeNXdy4FXo8tB+kdgZavlsNUH8DPgBXc/GZhApN7Q1GlmQ4HbgAp3H0vk8QZXhaDGx4CL2qxrt6bov82rgFOj7/lF9LN1dO4eui/gDODFVst3AXcFXVc7dT4LnA+sBgZH1w0GVgdYU0n0H8ZMYF50XWjqi9aQC6wnOtDdan1o6uTvz0fNI/IohHnABWGoESgFPurs99b2c0Pk2S1ndLb/ULYUiPOBtUEys1JgEvA3oNijT8SK/lkUXGX8FPg+0NxqXZjqAxgB1AK/inZzfmlm2YSoTnffAtwHbAI+IfLUs5fCVGMrHdV0XJ+jsIZCXA+sDYqZ9QOeBr7j7ruDrqeFmV0K1Lj7kqBr6UQKcBrwkLtPAvYRji7NEdF++eVAGTAEyDazrwdb1TE7rs9RWEMhtA+sNbNUIoHwO3f/U3T1djMbHH19MFATUHmfBy4zsw3Ak8BMM/ttiOprUQ1Uu/vfostPEQmJMNX5BWC9u9e6+2HgT8CZIauxRUc1HdfnKKyhEM9DbbudmRnwn8BKd/9Jq5eeA74R/f4bRMYaup273+XuJe5eSuR3tsDdvx6W+lq4+zZgs5mNjq46D1hBuOrcBEwzs6zo3/t5RAZDw1Rji45qeg64yszSzawMKAfe6XRvQQ3kxDGYcgmwBvgY+Oeg64nWdBaR5tcyYGn06xIgn8jg3tron3khqHUGfx9oDGN9E4HK6O/yGWBg2OoE/gVYBXwEPA6kB10j8ASRMY7DRFoCNx6tJuCfo5+h1cDF8RxDlzmLSIywdh9EJCAKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRj/H8M0zzBvxK6+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer at 0x7f6bff295198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAEICAYAAABWCOFPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAd4klEQVR4nO3deXhV9b3v8fc380ACZGQImFACiowaBq0oYh2x2l7trdZKHVquWo+n59ZW7Xlu5Xie02uf41NbT61eb4+1ta32VnvUInVEKjhUgyLKHBkDQgZkJkCS7/1j79DsmJAN7GStJJ/X8+Qha+211/omsD/8fr+11m+ZuyMi0iIp6AJEJFwUCiISQ6EgIjEUCiISQ6EgIjEUCiISQ6Eg3cbMrjGzl4KuQ47OdJ1C32RmDpS7e1UX7b8UWA+kuntjVxxDuoZaCnJczCw56BqkaygUegEzu8PMtpjZHjNbbWbnmdkUM3vLzHaa2Sdm9nMzS4tu/3r0rR+Y2V4z+6qZXWdmi9vs181sZPT7x8zsITObb2b7gHPNbJaZvW9mu81ss5nNbfX2lmPsjB7jjLbHMLMzzexdM9sV/fPMVq8tNLN/NbM3oj/XS2ZW0AW/PmlDodDDmdlo4FZgsrvnABcCG4Am4J+AAuAM4DzgFgB3Pzv69gnu3s/d/xDn4b4G/BuQAywG9gGzgQHALOBmM/tSdNuWYwyIHuOtNnXnAc8DDwD5wE+A580sv83xrgeKgDTg9jjrlBOgUOj5moB0YIyZpbr7Bnf/2N2XuPvb7t7o7huA/wOcc4LHetbd33D3ZndvcPeF7v5hdHkZ8MQxHGMWsNbdH4/W+ASwCvhiq21+5e5r3P0A8P+AiSdYv8RBodDDRQcKvwPMBWrM7EkzG2Jmo8xsnpltM7PdwI+ItBpOxObWC2Y21cxeM7NaM9sF3HQMxxgCbGyzbiMwtNXytlbf7wf6HWO9chwUCr2Au//e3c8CTgIc+DHwEJH/ecvdPRf4AWBH2c0+IKtlwcwGtXeoNsu/B54Dhrl7f+DhVsfo7LTW1mi9rQ0HtnTyPuliCoUezsxGm9lMM0sHGoADRLoUOcBuYK+ZnQzc3Oat24ERrZY/AE41s4lmlkGk5dGZHGCHuzeY2RQiYwAtaoHmNsdobT4wysy+ZmYpZvZVYAwwL47jShdSKPR86cC9QB2R5nYRkVbB7UQ+pHuA/wu0HUycC/w6enbiv7v7GuAe4BVgLZGBxM7cAtxjZnuAHxLp9wPg7vuJDEq+ET3GtNZvdPd64FLgu0A98H3gUnevi/9Hl66gi5dEJIZaCiISQ6EgIjEUCiISQ6EgIjFSgjpwQUGBl5aWBnV4kT5ryZIlde5e2NHrgYVCaWkplZWVQR1epM8ys7ZXksZQ90FEYigURCSGQkFEYgQ2ptCew4cPU11dTUNDQ9Cl9CoZGRmUlJSQmpoadCnSA4QqFKqrq8nJyaG0tBSzo93QJ/Fyd+rr66murqasrCzocqQHCFX3oaGhgfz8fAVCApkZ+fn5an1J3EIVCoACoQvodyrHInShICLBUiicgMcee4ytW7cmbH8bNmzg97///ZHlyspKbrvttoTtXyQeCoUTcDyh0NjY8XNR2oZCRUUFDzzwwHHXJ3I8FApt7Nu3j1mzZjFhwgTGjh3LH/7wB+655x4mT57M2LFjmTNnDu7OU089RWVlJddccw0TJ07kwIEDlJaWUlcXmTiosrKSGTNmADB37lzmzJnDBRdcwOzZs9mwYQPTp0/ntNNO47TTTuPNN98E4M4772TRokVMnDiR+++/n4ULF3LppZcCsGPHDr70pS8xfvx4pk2bxrJly47s+4YbbmDGjBmMGDFCISInLFSnJFv7lz8vZ8XW3Qnd55ghudz9xVOPus0LL7zAkCFDeP755wHYtWsX559/Pj/84Q8BuPbaa5k3bx5XXnklP//5z7nvvvuoqKjo9NhLlixh8eLFZGZmsn//fl5++WUyMjJYu3YtV199NZWVldx7773cd999zJsXmaZw4cKFR95/9913M2nSJJ555hkWLFjA7NmzWbp0KQCrVq3itddeY8+ePYwePZqbb75Z1yTIcVNLoY1x48bxyiuvcMcdd7Bo0SL69+/Pa6+9xtSpUxk3bhwLFixg+fLlx7zfyy67jMzMTCBykda3vvUtxo0bx1e+8hVWrFjR6fsXL17MtddeC8DMmTOpr69n165dAMyaNYv09HQKCgooKipi+/btx1yfSItOWwpm9iiRCTZr3H1sO68b8DPgEiJz81/n7u+daGGd/Y/eVUaNGsWSJUuYP38+d911FxdccAEPPvgglZWVDBs2jLlz53Z4zj8lJYXm5maAz2yTnZ195Pv777+f4uJiPvjgA5qbm8nIyOi0rvbm0mw51Zienn5kXXJy8lHHLUQ6E09L4THgoqO8fjFQHv2aQ+R5Az3W1q1bycrK4utf/zq33347770XybeCggL27t3LU089dWTbnJwc9uzZc2S5tLSUJUuWAPD00093eIxdu3YxePBgkpKSePzxx2lqamp3f62dffbZ/O53vwMi3YqCggJyc3NP7IcVaUenLQV3fz36WPGOXA78xiP/lb1tZgPMbLC7f5KgGrvVhx9+yPe+9z2SkpJITU3loYce4plnnmHcuHGUlpYyefLkI9ted9113HTTTWRmZvLWW29x9913c+ONN/KjH/2IqVOndniMW265hSuuuII//vGPnHvuuUdaEePHjyclJYUJEyZw3XXXMWnSpCPvmTt3Ltdffz3jx48nKyuLX//61133S5A+La4p3qOhMK+D7sM84F53XxxdfhW4w90/M4OKmc0h0ppg+PDhp2/cGDvXw8qVKznllFOO/aeQTul3Ky3MbIm7dzg6noiBxvauoW03adz9EXevcPeKwsIOZ4MSkQAlIhSqgWGtlkuIPCdQRHqgRITCc8Bsi5gG7DqR8QQ9sSrx9DuVYxHPKckngBlAgZlVA3cDqQDu/jCRB4VeAlQROSV5/fEWk5GRQX19vW6fTqCW+RTiOe0pAvGdfbi6k9cd+HYiiikpKaG6upra2tpE7E6iWmZeEolHqC5zTk1N1exAIgHTZc4iEkOhICIxFAoiEiPUobB62x5+8tJqnVIT6UahDoVl1Tt5YEEVH21J7LwKItKxUIfCF04pJjnJeHH5tqBLEekzQh0KA7PTmFKap1AQ6UahDgWAC04tZm3NXtbV7g26FJE+oQeEwiAAXlyuKcZEukPoQ2HogEzGDe3PSyvUhRDpDqEPBYALTy3m/U07WbO9/anKRCRxekQofHXycPKz07jtifdpONwUdDkivVqPCIXCnHTu+8oEVm3bw71/WRV0OSK9Wo8IBYBzTy7ihs+X8dibG/jOk+9Ts1uPVhfpCqG6dbozd11yMv3Sk3n4r+t4ZWUNN5xVxo2fL6N/lp6GJJIocc3m3BUqKiq8svIzEz7HZX3dPn78l1W8sHwbOekpXDPtJK47s5RB/TW7kEhnOpvNuUeGQosVW3fz4GtV/OWjT0hOMi4aO5jZZ5xExUkDNZ2bSAd6dSi02LxjP796YwN/XLKZPQ2NjC7O4aopw/jypKEMyEpLyDFEeos+EQot9h9q5NmlW3ninU0sq95FWnIS548p5orThzK9vJDU5B4zrirSZfpUKLS2fOsu/lhZzbNLt/Dp/sPkZ6fxxQlDuHziECYOG6DuhfRZfTYUWhxqbGbh6hr+6/0tvLqqhkONzQzPy+LS8YO5dPwQThmco4CQPqXPh0Jruw4c5sWPtvHnZVt58+N6mpqdsoJsZo0bzMXjBjFmcK4CQno9hUIH6vce5IXl23h+2Se8va6eZofheVlcNHYQF55azKRhA0lKUkBI76NQiEP93oO8tGI7Ly7fxhtVdRxucgr6pXP+mCLOH1PMmZ8rICM1OegyRRJCoXCMdjcc5rVVNby0fDsLV9ew71ATWWnJnDWygC+MKebc0UUU5qQHXabIcessFHrUZc7dITcjlcsnDuXyiUM52NjEWx/X88rK7by6soaXVkQmepkwbAAzRxcx8+QiTh2Sq26G9CpqKcTJ3VnxyW4WrKzh1VU1fFC9E/fIHZwzRhUyY3QRZ5UX0D9T92FIuCWk+2BmFwE/A5KBX7r7vW1e7w/8FhhOpPVxn7v/6mj77Gmh0Fbd3oMsXF3LwtU1vL6mlt0NjSQnGacNH8A5owo5Z5RaERJOJxwKZpYMrAHOB6qBd4Gr3X1Fq21+APR39zvMrBBYDQxy90Md7benh0JrjU3NLN28MxISa2qOPKciPzuN6eUFTC8vZPqoAopydMOWBC8RYwpTgCp3Xxfd4ZPA5cCKVts4kGORk/z9gB1A43FX3cOkJCdRUZpHRWket184mto9B1lcVctfV9eyaG0dzyzdCsApg3M5OxoSFaUDdUZDQimelsKVwEXu/s3o8rXAVHe/tdU2OcBzwMlADvBVd3++nX3NAeYADB8+/PSNGzcm6ucIreZmZ/nW3SyqquX1NbUs2fgph5uc9JQkppTlMb28gLNGFurKSuk2iWgptPcvtW2SXAgsBWYCnwNeNrNF7h7zvDd3fwR4BCLdhziO3eMlJRnjSvozrqQ/t8wYyb6DjfxtfT2L1taxaG0dP5q/ClhFQb80Pj+ygLNGRloSmhtCghJPKFQDw1otlwBb22xzPXCvR5odVWa2nkir4Z2EVNmLZKenMPPkYmaeXAzAJ7sOsHhtHYur6nijqo5no12NkUX9OGtkAZ8fWcC0EXnkZOishnSPeLoPKUQGGs8DthAZaPyauy9vtc1DwHZ3n2tmxcB7wAR3r+tov71poDFRmpudVdv2sLiqlsVV9byzvp6Gw82kJBkThw3grPJIS2LCsAG6DVyOW6JOSV4C/JTIKclH3f3fzOwmAHd/2MyGAI8Bg4l0N+51998ebZ8Khc4dbGxiycZPj7QkPtyyC3fol57CtBH5nD0q0tUozc/SeITETZc59yI79x/izY8j4xGLq2rZvOMAACUDM5leXsjZ5QWcOVIXUMnRKRR6sQ11+1hUVcfra2p56+N69h6MXEA1aVjkAqoZo3UBlXyWQqGPONzUzPubdvL6mlpeX1vLsupdABT0S2fG6ELOHV3E9FEF5GrAss9TKPRRdXsP8vqaWl5bHbk+YteBw6QkGVNH5PGFU4o5f0wxJQOzgi5TAqBQEBqbmnlv004WrKrhlZXbqarZC8DYoblcOGYQl4wfzOcK+wVcpXQXhYJ8xvq6fby0fBsvLt/Ge5t2AnDyoBwumziEyyYMUQuil1MoyFFt29XAXz76hD9/sPVIQEwty+PK00uYNX4wWWmacqO3UShI3Dbv2M8z72/h6feq2VC/n5z0FL582lBmn3ESI4tygi5PEkShIMfM3Xl3w6c88c4mnv/wEw41NnPu6EJuOudzTB2RH3R5coIUCnJC6vYe5Hdvb+LxtzdQt/cQU8sit4dPLs0LujQ5TgoFSYgDh5p44p1NPPzXj6nZc5CLxw7iB5ecwrA8DUr2NJ2Fgu6qkbhkpiVzw1ll/PV75/I/zx/FwtW1XHD/6zy6eD1NzX3iLvg+Q6EgxyQzLZnbzivn1e+ew9QRedwzbwXX/uffqN97MOjSJEEUCnJchgzI5FfXTebHV4yjcuOnfPE/FrOsemfQZUkCKBTkuJkZX508nD/dfCZmxtWPvM27G3YEXZacIIWCnLCxQ/vz9M1nUtw/g288+g5vr6sPuiQ5AQoFSYhB/TN4cs40hgzIZM5vKtlUvz/okuQ4KRQkYYpyMnj0G5MxM/7Hb5dw4FBT0CXJcVAoSEINz8/ip1dNZNW23cx9bnnnb5DQUShIwp07uog5Z4/gD5WbeW/Tp0GXI8dIoSBd4h9mllOYk849f15Bsy5u6lEUCtIl+qWn8P0LR7N0806e/WBL0OXIMVAoSJe54rQSxg3tz09eXqPWQg+iUJAuk5RkfHN6GZt3HOCNjzt8LpCEjEJButSFpw5iQFYqT76zOehSJE4KBelSGanJ/LdJJby0YptumuohFArS5a6eMozDTc7T71UHXYrEQaEgXa68OIfTTxrIU0sUCj1BXKFgZheZ2WozqzKzOzvYZoaZLTWz5Wb218SWKT3dhacWs2b7Xmp2NwRdinSi01Aws2TgQeBiYAxwtZmNabPNAOAXwGXufirwlS6oVXqwadEJX9/SHZShF09LYQpQ5e7r3P0Q8CRweZttvgb8yd03Abh7TWLLlJ5uzOBcctJTeHud5lsIu3hCYSjQ+nxSdXRda6OAgWa20MyWmNns9nZkZnPMrNLMKmtra4+vYumRUpKTmFKWx9/UUgi9eEKhveeYt708LQU4HZgFXAj8LzMb9Zk3uT/i7hXuXlFYWHjMxUrPNm1EPuvq9rFd4wqhFk8oVAPDWi2XAFvb2eYFd9/n7nXA68CExJQovUXLuIJmZgq3eELhXaDczMrMLA24CniuzTbPAtPNLMXMsoCpwMrElio93ZghueRkpCgUQq7Tp4e6e6OZ3Qq8CCQDj7r7cjO7Kfr6w+6+0sxeAJYBzcAv3f2jrixcep7kJGNqWZ4GG0MurkcKu/t8YH6bdQ+3Wf534N8TV5r0RmOH9ufVVTU0HG4iIzU56HKkHbqiUbpVWUE27rBRE7uGlkJButWIgn4ArK/bF3Al0hGFgnSr0oLIA2kVCuGlUJBulZORSkG/dNbX7Q26FOmAQkG63YiCbLUUQkyhIN2urCCb9XUaaAwrhYJ0u9KCbOr2HmR3w+GgS5F2KBSk25UVZAOwQV2IUFIoSLcbURgJBY0rhJNCQbrd8LwszBQKYaVQkG6XkZrMkP6ZCoWQUihIIEYU6rRkWCkUJBCl+QqFsFIoSCAG9c9gT0MjDYebgi5F2lAoSCDystMA2LHvUMCVSFsKBQmEQiG8FAoSiJZQqFcohI5CQQLREgqfKhRCR6EggchXSyG0FAoSiNyMVJKTjB379Hj6sFEoSCCSkoyBWans2Kc7JcNGoSCByctOU0shhBQKEphIKGhMIWwUChIYhUI4KRQkMAqFcFIoSGDystPZeeAwTc1tH2IuQVIoSGDyslJxh5371VoIk7hCwcwuMrPVZlZlZnceZbvJZtZkZlcmrkTprfL6pQO6/yFsOg0FM0sGHgQuBsYAV5vZmA62+zGRp1OLdEpXNYZTPC2FKUCVu69z90PAk8Dl7Wz3D8DTQE0C65NebGCW7n8Io3hCYSiwudVydXTdEWY2FPgyEPN4+rbMbI6ZVZpZZW1t7bHWKr1Mfj+1FMIonlCwdta1HS7+KXCHux91Gh13f8TdK9y9orCwMN4apZdqaSloTCFcUuLYphoY1mq5BNjaZpsK4EkzAygALjGzRnd/JiFVSq+UlpJETkaKQiFk4gmFd4FyMysDtgBXAV9rvYG7l7V8b2aPAfMUCBIPXcAUPp2Ggrs3mtmtRM4qJAOPuvtyM7sp+vpRxxFEjkahED7xtBRw9/nA/Dbr2g0Dd7/uxMuSviI/O40tOxuCLkNa0RWNEqiBWWk6JRkyCgUJVF6/SPfBXfc/hIVCQQI1MCuNQ03N7D+kh8KEhUJBApWbkQrAnobGgCuRFgoFCVRuZmSse3eD5moMC4WCBCrnSEtBoRAWCgUJVG5GtKVwQN2HsFAoSKBaWgrqPoSHQkEC9fcxBbUUwkKhIIFqOfuw+4BaCmGhUJBAZaQmk5acpO5DiCgUJHC5mSm6TiFEFAoSuNyMVHUfQkShIIHLyVBLIUwUChK43MxUjSmEiEJBAqfuQ7goFCRw6j6Ei0JBAqfuQ7goFCRwuRkpNBxu5lBjc9ClCAoFCQHdKRkuCgUJnO5/CBeFggQuJ133P4SJQkECl5upKdnCRKEggdOUbOGiUJDA5ej26VBRKEjgWqZkU/chHBQKErjstBSSTN2HsIgrFMzsIjNbbWZVZnZnO69fY2bLol9vmtmExJcqvVVSktEvPUXdh5DoNBTMLBl4ELgYGANcbWZj2my2HjjH3ccD/wo8kuhCpXfLzUxV9yEk4mkpTAGq3H2dux8CngQub72Bu7/p7p9GF98GShJbpvR2uRm6/yEs4gmFocDmVsvV0XUduRH4S3svmNkcM6s0s8ra2tr4q5ReLycjRc9+CIl4QsHaWdfuI4LN7FwioXBHe6+7+yPuXuHuFYWFhfFXKb2e7pQMj5Q4tqkGhrVaLgG2tt3IzMYDvwQudvf6xJQnfUVuhsYUwiKelsK7QLmZlZlZGnAV8FzrDcxsOPAn4Fp3X5P4MqW3i3Qf1FIIg05bCu7eaGa3Ai8CycCj7r7czG6Kvv4w8EMgH/iFmQE0untF15UtvU1uZip7DzXS3OwkJbXXY5XuEk/3AXefD8xvs+7hVt9/E/hmYkuTvqR/ZirusOvAYQZmpwVdTp+mKxolFErzswBYV7cv4EpEoSChUF6UA8Da7XsCrkQUChIKJQMzyUhNYm3N3qBL6fMUChIKSUnGyKJ+rFFLIXAKBQmNUUU5rN2ulkLQFAoSGuXFOWzb3aArGwOmUJDQKC/qB6DWQsAUChIao4ojZyCqajSuECSFgoRGyxmINWopBEqhIKHRcgZCpyWDpVCQUCkvytEFTAGL694Hke5SXtyP/3p/Cz95eQ3TywvITE0OuqQe5dQhuURvSjxuCgUJlYvHDmbByhr+Y8FaHnh1bdDl9Djr//clJ7wPhYKESllBNk/dfCY79h1i6eZPaWxqd5Iv6UIKBQmlvOw0Zp5cHHQZfZIGGkUkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkhkJBRGIoFEQkRlyhYGYXmdlqM6syszvbed3M7IHo68vM7LTElyoi3aHTUDCzZOBB4GJgDHC1mY1ps9nFQHn0aw7wUILrFJFuEk9LYQpQ5e7r3P0Q8CRweZttLgd+4xFvAwPMbHCCaxWRbhBPKAwFNrdaro6uO9ZtMLM5ZlZpZpW1tbXHWquIdIN4QqG9WSDbzpEVzza4+yPuXuHuFYWFhfHUJyLdLJ5QqAaGtVouAbYexzYi0gPEEwrvAuVmVmZmacBVwHNttnkOmB09CzEN2OXunyS4VhHpBp1O3OrujWZ2K/AikAw86u7Lzeym6OsPA/OBS4AqYD9wfdeVLCJdKa7ZnN19PpEPfut1D7f63oFvJ7Y0EQmCrmgUkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRgKBRGJYZF7mQI4sFktsLGTzQqAum4o50SoxsRQjYkRT40nuXuHsxwFFgrxMLNKd68Iuo6jUY2JoRoTIxE1qvsgIjEUCiISI+yh8EjQBcRBNSaGakyME64x1GMKItL9wt5SEJFuplAQkRihDYXOHmobBDMbZmavmdlKM1tuZv8YXZ9nZi+b2dronwMDrjPZzN43s3lhrC9a0wAze8rMVkV/n2eErU4z+6fo3/NHZvaEmWUEXaOZPWpmNWb2Uat1HdZkZndFP0OrzezCeI4RylCI86G2QWgEvuvupwDTgG9H67oTeNXdy4FXo8tB+kdgZavlsNUH8DPgBXc/GZhApN7Q1GlmQ4HbgAp3H0vk8QZXhaDGx4CL2qxrt6bov82rgFOj7/lF9LN1dO4eui/gDODFVst3AXcFXVc7dT4LnA+sBgZH1w0GVgdYU0n0H8ZMYF50XWjqi9aQC6wnOtDdan1o6uTvz0fNI/IohHnABWGoESgFPurs99b2c0Pk2S1ndLb/ULYUiPOBtUEys1JgEvA3oNijT8SK/lkUXGX8FPg+0NxqXZjqAxgB1AK/inZzfmlm2YSoTnffAtwHbAI+IfLUs5fCVGMrHdV0XJ+jsIZCXA+sDYqZ9QOeBr7j7ruDrqeFmV0K1Lj7kqBr6UQKcBrwkLtPAvYRji7NEdF++eVAGTAEyDazrwdb1TE7rs9RWEMhtA+sNbNUIoHwO3f/U3T1djMbHH19MFATUHmfBy4zsw3Ak8BMM/ttiOprUQ1Uu/vfostPEQmJMNX5BWC9u9e6+2HgT8CZIauxRUc1HdfnKKyhEM9DbbudmRnwn8BKd/9Jq5eeA74R/f4bRMYaup273+XuJe5eSuR3tsDdvx6W+lq4+zZgs5mNjq46D1hBuOrcBEwzs6zo3/t5RAZDw1Rji45qeg64yszSzawMKAfe6XRvQQ3kxDGYcgmwBvgY+Oeg64nWdBaR5tcyYGn06xIgn8jg3tron3khqHUGfx9oDGN9E4HK6O/yGWBg2OoE/gVYBXwEPA6kB10j8ASRMY7DRFoCNx6tJuCfo5+h1cDF8RxDlzmLSIywdh9EJCAKBRGJoVAQkRgKBRGJoVAQkRgKBRGJoVAQkRj/H8M0zzBvxK6+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u\n",
    "*dfw(sw.faceValue)*[[1.0]]) + u*(fw(sw.faceValue)-dfw(sw.faceValue)*sw.faceValue).divergence == 0\n",
    "\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "steps = 50\n",
    "viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)\n",
    "for step in range(steps):\n",
    "    sw.updateOld()\n",
    "    swres = 1.0e6\n",
    "    while swres > 1e-5:\n",
    "        swres = eq.sweep(dt = dt, var = sw)\n",
    "#         print(swres)\n",
    "        \n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fractional_flow as ff\n",
    "xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil = ff.frac_flow_wf(muw=muw, muo=muo, ut=u, phi=1.0, \\\n",
    "  k=1e-12, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0, \\\n",
    "  nw=nw, sw0=swc, sw_inj=1.0, L=L, pv_inj=5.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcp0lEQVR4nO3deXRc5Z3m8e9PuyxZm1WSbcmybGMbS8YGLAwYGzDpsLhDnHQ2IN10usMQzmk6PWcmfUJPJpk5J9N9kk736U4mJB46ISQhgU5DOiFgFp8EzGqwjPG+yfImy4u8SV4ka3vnj1sSsqylSirVvVV6PufUKdW9b937eynq8av3LjLnHCIikvhS/C5ARERiQ4EuIpIkFOgiIklCgS4ikiQU6CIiSSLNrx0XFxe7yspKv3YvIpKQNmzYcMI5FxponW+BXllZSW1trV+7FxFJSGZ2YLB1mnIREUkSCnQRkSShQBcRSRIKdBGRJKFAFxFJEsMGupk9bmbHzWzrIOvNzL5nZnVmttnMro19mSIiMpxIRuhPAHcOsf4uYHb48SDww9GXJSIi0Ro20J1zrwOnhmiyEviZ86wDCsxsSqwK7K+to4tnNjSg2/6KiFwqFnPoZcChPq8bwssuY2YPmlmtmdU2NTWNaGffenEnX/mPTazdPbL3i4gkq1gEug2wbMDhs3PuMedcjXOuJhQa8MrVYTWdvQjA2bbOEb1fRCRZxSLQG4BpfV6XA40x2O7Awv98aMJFRORSsQj054D7w2e73AA0O+eOxGC7A+r5dUBz6CIilxr25lxm9hRwK1BsZg3A/wLSAZxzq4DVwAqgDrgA/MVYFRuuB2/fY7kXEZHEM2ygO+fuHWa9A/4qZhUNI6V3ykWJLiLSV8JdKfrhlIuvZYiIBE7iBbqmXEREBpR4gR5+Vp6LiFwq4QK997RFDdFFRC6RcIFu4URXnIuIXCrxAl1zLiIiA0q8QAcyaaei4Tlo/MDvckREAiPxAj08Ql+09f/A+h/5W4yISIAkXqBjXCSDw6XLYcfvoLPd75JERAIh8QI9PELfP+UuaDsD9a/6W5CISEAkbKAfKb4RsvJh66/9LUhEJCASLtB7TkTvsnSYdzfsfAE62nyuSUTEfwkX6Nb3wqLqP4H2s1C3xt+iREQCIPECPfzsAGbcAhMmadpFRIREDPTeETqQmgZVK2H3S3DxrK91iYj4LeECPSM1Ffjwb4ty9Z9CRyu88BXdglFExrWEC/T7b5xOVnoK33+1jru++wY/3l/E+SV/C5uf1oVGIjKuJVygVxbnsO7vPsI3V1aTkWp88/ntXP3aQjZl30D3i4/Qse8dv0sUEfGF+XUb2pqaGldbWzvq7ew+dpZnNzTwyvu7+Un731Kc0sLbU+6nYPnfUHPFVFJ6/madiEgSMLMNzrmaAdcleqD36OzqZv3GjeS8+j9ZcP5tDrtJPJH2WWzhPay4ppKF5fm9f+1IRCRRjYtA76t192u0rf4ahWe2cswV8njnnbyRdxe3LJzL3QumMm/KRIW7iCSkcRfogHfGS/1rdLz+L6QfWEu7ZfCbziX8vPOPaC6o5o75k7mjejLXVhRqWkZEEsb4DPS+jm6F9f+G2/QrrPMCB9Jn8tO2ZTzbcRPpuZP4aFUpt1eXsmTWJDLTUuNTk4jICCjQe7SegS3/ARt/Dkc20ZWSzubsG3j87PW81L6ArMwsbr2yhNurSrl5Toj87PT41iciMgwF+kCObIJNT3sBf76Jjox8NuYs48fNi1hzYTaWksqi6YUsn1vCbVeWMKc0V/PuIuI7BfpQujq9e6pveQZ2Pg/t5+jIKmZr/i08db6GZ05Mo5sUygqyWX5liOVzS1gyq5jsDE3NiEj8KdAj1X4B9rwM2/4Tdr8Cna10TSimftJynu+o4YnGcprbjcy0FG6cNYlb5oRYNruYWSGN3kUkPhToI3HxnHdb3u2/9cK94zwuM4+mKbeyNmUxjx+dxY5T3n+7KflZLL2imKWzi1l6RTGTcjN9Ll5EkpUCfbQ6WmHvq96UzK7V0HoaUjNoLb+JrTlL+M2Fq3j+QCrNrR0AVE3JY9nsYpbNDlFTWUhWuqZnRCQ2FOix1NUJh9bBztVeuJ/eB4Arnc+xybfypl3LM0dL2XCohY4uR2ZaCotnFLFkVjE3zCziqrJ80lIT7hY6IhIQCvSx4hyc2A27XoQ9r8DBdeC6ILuIzpm3sWvijbzQOo81+zrZc/wcADkZqdRUFnHjrEncMHMS86fmKeBFJGIK9HhpPQ11v4c9a7z59wsnAYOyRZyvuJWNmYt4+fRU3tnXTF044HMz06ipLOTGmV7AVyvgRWQICnQ/dHdB4wdesO9ZA4c3AA6yCmDmrbSULePdlGt47WgG6+pPsrfpPAATM9O4bkYRi2cUcV1lIfPL8nX1qoj0UqAHwYVT3vnudb+HvX+As0e85ZNmw6zbODPlJt7qquKtQ22s23uS+hNewGempbBwWgHXVRZSU1nEtRWFuoJVZBxToAeNc3B8hxfs9a/C/regsxVS0qCsBmbeyunJS3ivYwbvHTxH7YHTbDvcTGe3wwzmlk7kusoiaioLua6yiKkF2X73SETiZNSBbmZ3At8FUoEfOee+1W99PvAkUAGkAf/knPvJUNsc14HeX0cbNLznnRpZ/6o3VYOD9ByYvgRm3kLrtKVsvFhG7YFm1u8/xfsHTnO+vQuAsoJsasIj+EUVhcwpzdU8vEiSGlWgm1kqsBv4KNAArAfudc5t79PmfwD5zrmvmlkI2AVMds61D7ZdBfoQWk/D/jeh/jWoXwsn93jLs4tgxjKoXEbn9GXs7JhM7YHTrD9wmvX7TnE8/IezJ2SksqA8n2srCrmmopBrKgoo1sVOIklhqEBPi+D9i4E651x9eGNPAyuB7X3aOGCiede/5wKngM5RVT2eZRfCvLu9B0BLI+x73XvUr4XtvyUNmJ87mfmVS/nC3JtxdyyjgVLeP3SGjQfPsPHgaR57vZ7Obu8f7GlF2V7ATyvgmopC5k3JIyNNo3iRZBLJCP3TwJ3OuQfCr/8MuN4593CfNhOB54ArgYnA55xzLwywrQeBBwEqKioWHThwIFb9GD+cg1P1sP+NcMi/AeePe+vyyqFyaXgUv5S2nHK2NrZ4AX/oNO8fOMPRljbAO9h6VVk+11R4AX9tRSGT87N87JiIRGK0Uy6fAe7oF+iLnXN/3afNp4GbgP8GzALWAAudcy2DbVdTLjHSc3HT/je8cN//Jlw44a3Ln+YFfM+jYDpHWtp6R/DvHzzDlsPNtHd2A949aRaWF7BgWj4Lywu4qjyfvCydUSMSJKOdcmkApvV5XQ409mvzF8C3nPevQ52Z7cMbrb83gnolGmYQmus9rnvAC/imnV6w73vdu4J101Ne27xyplQuZUrlTay4cSmsuJH2LseOIy29Ab+p4QwvbTvau/mZxTksnFbAgvJ8FpQXUD01T/emEQmoSEboaXgHRT8CHMY7KHqfc25bnzY/BI455/63mZUC7+ON0E8Mtl2N0OOkb8DvfxMOvAXnm7x1E6dC5U3hEfwyKJoJZpy50M7mhmY2N5xhU/j5WIt3wDUtxZhTOpGF07yAX1Cez5zSiaTrrBqRuIjFaYsrgH/FO23xcefc35vZQwDOuVVmNhV4ApgCGN5o/cmhtqlA90nfKZr9b3rnwPfMwedO/jDgpy+F4tnebwDA0eY2NjWcYXPDmXDYN/feXTIzLYXqqXksKC/oDfoZk3L0x7dFxoAuLJLBOQcn6z4cwe9/E86Fp1xySvrMwS+7JOCdcxw4eSEc8s1sOnSGrY3NtHV48/ETs9KYPzWfq8rzqZ6ax1Vl+VQq5EVGTYEuket7Fk3PCP5s+JDJEAEP0NnVzZ7j53qnarYebmbnkbO0d3khn5uZRlU43OeXec8zinNJVciLREyBLiPXE/AH3gqfRfPGh/ehyZ0cPk3yZu9RWHlJwAN0dHWz+9hZth1uYcvhZrY2NrO9sYWL4TNrJmSkUjUlj/ll+cwvy+eqsnxmhXJ0pavIIBToEjuXnAcfDvhzx7x1+RUw82aYcYv3mFg64CY6u7qpazrH1sMtbD3sjeS3NbbQ2uHdyiArPYV5U/K8KZuyfKrL8nTgVSRMgS5jp+cg677XYd9aL+TbznjrQvNg5q0wazlMvwkycwfdTFe3Y9+Jc94oPjya397YwrmL3gXHGWkpzJs8kerwKH7+1HzmTM7VrYVl3FGgS/x0d8HRzd4tCupfg4PvQGcbpKRDxQ0w6za44iMwecFl0zOXbarbsf/kebaER/BbGrwpm7NtXsinpxpXlEykemoe86fmUV2Wz7wpeeRmRnJ5hUhiUqCLfzravFDf+wfvbpLHtnjLc0vhij+C2bd7IZ+VF9HmnHMcPHWhN+S3NbawvbGZE+e8+8CZQeWkHKqm5lE9NY/qqd5ZNro5mSQLBboEx9mjXrjvWeM9t53x7gM/fQnMXeE9CqdHtUnnHMfPXuydi9/W6D03nG7tbTM5Lysc8HlUhUO+vDAbG+a3BJGgUaBLMHV1eveB3/2y94e2T+zylpdeBVUroerj3i0NRqj5Qgfbjnhz8dsavQOwe5vOEb4BJfnZ6eEzbD4cyc8M6TRKCTYFuiSGk3th12rY8Ts49K63rKQa5v8JXPVp77TIUWpt72Ln0ZZLpmt2HD3be4OyrPQU5k/N558/u5Dpk3JGvT+RWFOgS+JpOQI7noOtz34Y7hVL4Or7oPoTkDkxZrvq7Opmb9N5th5u5t19J/lVbQPfvedqVl5dFrN9iMTKUIGuE3slmPKmwPVfgi++An+zGW77unfPmecehn++En73X+HolpjsKi01hbmTJ/KpReU8dMusmGxTxA8KdAm+wulw81fg4Vr44hqo+oR3S+BVS+EnK7z59+5uv6sU8Z0CXRKHGUxbDJ94FP77Trj97+HMQXjqHvjBDbDlGQW7jGsKdElM2YWw5GH48kb41I8hJRWe/SL8v2XeWTMi45ACXRJbarp3BsxDb3nB3tEKv/ws/PJzcGqf39WJxJUCXZJDSooX7H/1rjcVs/9Nbxrmre9qGkbGDQW6JJfUdG8q5uH1MPujsOYb8OQnvStURZKcAl2SU95U+OzP4e7vwcF34Yc3wcF1flclMqYU6JK8zGDRn8OX1kJ2AfxsJexc7XdVImNGgS7JLzQX/vJlKK2Gf/88vP8zvysSGRMKdBkfcorh/udg5nJ47sve/WJEkowCXcaPzFy45xdQtgie/S/QuNHvikRiSoEu40t6Ntz7FOSE4Jf3QHOD3xWJxIwCXcaf3BL4/K+g4wI8+4DOU5ekoUCX8alkHtzxD96fx/vgSb+rEYkJBbqMX9f8qXeP9TXfgPMn/K5GZNQU6DJ+mcHH/gUunoNXvu53NSKjpkCX8a3kSrjpy7Dpl94VpSIJTIEusuwr3u143/m/flciMioKdJGMCXDt/bDzBVLPHva7GpERU6CLAFz3AAB5W3VbAElcCnQRgIIKmLuCvO2/IJN2v6sRGREFukiP679EattpPp76tt+ViIyIAl2kR+Uy2ouu5Aup+pukkpgU6CI9zGipuo/qlANMOH/Q72pEohZRoJvZnWa2y8zqzOyRQdrcamYfmNk2M1sb2zJF4uPC9OUAlBzXtIsknmED3cxSgUeBu4Aq4F4zq+rXpgD4AfBx51w18JkxqFVkzHXmz+BQd4iS42/5XYpI1CIZoS8G6pxz9c65duBpYGW/NvcBv3bOHQRwzh2PbZkicWLGG93zKW56F7o6/a5GJCqRBHoZcKjP64bwsr7mAIVm9pqZbTCz+wfakJk9aGa1Zlbb1NQ0sopFxtgb3QtI7zwHhzf4XYpIVCIJdBtgmev3Og1YBPwxcAfwdTObc9mbnHvMOVfjnKsJhUJRFysSD293V+Mw2PsHv0sRiUokgd4ATOvzuhxoHKDNS8658865E8DrwMLYlCgSX83kcqZwPtS/6ncpIlGJJNDXA7PNbIaZZQD3AM/1a/NbYJmZpZnZBOB6YEdsSxWJn+Ohm6ChFtqa/S5FJGLDBrpzrhN4GHgZL6R/5ZzbZmYPmdlD4TY7gJeAzcB7wI+cc1vHrmyRsdVUsgRcF+x7w+9SRCKWFkkj59xqYHW/Zav6vf4O8J3YlSbin1NFCyAj15tHn/cxv8sRiYiuFBUZgEvJgLJF0Pi+36WIREyBLjKY0Fw4sQdc/5O6RIJJgS4ymOI50H4OWvqf1CUSTAp0kcGE5nrPJ3b5W4dIhBToIoMpDgd6025/6xCJkAJdZDC5JZCVrxG6JAwFushgzLxR+ok9flciEhEFushQiudAk0bokhgU6CJDCc2B88eh9bTflYgMS4EuMhQdGJUEokAXGUoofBfoEwp0CT4FushQCqZDaqbOdJGEoEAXGUpKKhTP1pSLJAQFushwimdrhC4JQYEuMpziuXD6AHS0+l2JyJAU6CLDCc0BHJzc63clIkNSoIsMp1g36ZLEoEAXGU5Bhfes2+hKwCnQRYaTOdE7dfF8k9+ViAxJgS4yHDPICcE5BboEmwJdJBK5IY3QJfAU6CKRyAl5N+kSCTAFukgkckJw/oTfVYgMSYEuEomc8JSLc35XIjIoBbpIJHJC0NUObc1+VyIyKAW6SCRyS7xnTbtIgCnQRSKRU+w968CoBJgCXSQSOSHvWacuSoAp0EUikdMz5aJAl+BSoItEYsIk71lXi0qAKdBFIpGaBtlFGqFLoCnQRSKVW6KDohJoCnSRSOlqUQk4BbpIpHJ0gy4JNgW6SKR0C10JuIgC3czuNLNdZlZnZo8M0e46M+sys0/HrkSRgMgJwcVm6LzodyUiAxo20M0sFXgUuAuoAu41s6pB2n0beDnWRYoEQq4uLpJgi2SEvhioc87VO+fagaeBlQO0+2vgWUCnAUhy0tWiEnCRBHoZcKjP64bwsl5mVgZ8Elg11IbM7EEzqzWz2qYmfSkkwfQGus50kWCKJNBtgGX9bwr9r8BXnXNdQ23IOfeYc67GOVcTCoUirVEkGHoC/Zx+CZVgSougTQMwrc/rcqCxX5sa4GkzAygGVphZp3PuNzGpUiQINOUiARdJoK8HZpvZDOAwcA9wX98GzrkZPT+b2RPA8wpzSTqZuZA+QYEugTVsoDvnOs3sYbyzV1KBx51z28zsofD6IefNRZJKTrECXQIrkhE6zrnVwOp+ywYMcufcF0ZflkhA6WpRCTBdKSoSjZwSXS0qgaVAF4mGplwkwBToItHICcGFE9Dd7XclIpdRoItEI7cEujuh7YzflYhcRoEuEo3sQu+59bS/dYgMQIEuEo2sAu9ZI3QJIAW6SDSy8r3nVgW6BI8CXSQa2RqhS3Ap0EWi0TPlohG6BJACXSQaGqFLgCnQRaKRng2pmRqhSyAp0EWilV2gEboEkgJdJFpZBRqhSyAp0EWipRG6BJQCXSRaGqFLQCnQRaKlEboElAJdJFpZBdDa7HcVIpdRoItEK7sALjZDd5fflYhcQoEuEq3eG3RplC7BokAXiVa2Al2CSYEuEi3dQlcCSoEuEq1s3aBLgkmBLhItjdAloBToItHSCF0CSoEuEi2N0CWgFOgi0UrPhpR0jdAlcBToItEy0+X/EkgKdJGR0A26JIAU6CIjoRG6BJACXWQkNEKXAFKgi4yERugSQAp0kZHQCF0CSIEuMhLZBd7Nubq7/a5EpJcCXWQksgoAB+1n/a5EpJcCXWQkdPm/BFBEgW5md5rZLjOrM7NHBlj/eTPbHH68bWYLY1+qSIDo8n8JoGED3cxSgUeBu4Aq4F4zq+rXbB9wi3NuAfBN4LFYFyoSKBqhSwBFMkJfDNQ55+qdc+3A08DKvg2cc287506HX64DymNbpkjAaIQuARRJoJcBh/q8bggvG8wXgRcHWmFmD5pZrZnVNjU1RV6lSNBohC4BFEmg2wDL3IANzZbjBfpXB1rvnHvMOVfjnKsJhUKRVykSNBqhSwClRdCmAZjW53U50Ni/kZktAH4E3OWcOxmb8kQCKiMHUtI0QpdAiWSEvh6YbWYzzCwDuAd4rm8DM6sAfg38mXNud+zLFAkYM2+UrhG6BMiwI3TnXKeZPQy8DKQCjzvntpnZQ+H1q4BvAJOAH5gZQKdzrmbsyhYJgGxd/i/BEsmUC8651cDqfstW9fn5AeCB2JYmEnBZ+RqhS6DoSlGRkcopgZYjflch0kuBLjJSpVVwcg90XvS7EhFAgS4ycqXV0N0JJ3QegASDAl1kpErne89Ht/pbh0iYAl1kpIpmQWomHFOgSzAo0EVGKjUNSq6EY9v8rkQEUKCLjE7pfAW6BIYCXWQ0Sqvh/HE4d9zvSkQU6CKj0nNgVKN0CQAFusholFZ7zwp0CQAFusho5BRD7mQFugSCAl1ktEqr4dgWv6sQUaCLjFppNTTtgq4OvyuRcU6BLjJapfOhqx1O1vldiYxzCnSR0eo5MLrxSWhu8LcWGdciuh+6yHjzD6t38P0/RDbiTnMdrEqZxvR3vg/vfJ8jVkK7ZfauH/AP8ErcrE1fxpOZ9/hdxiU+d900Hlg2M+bbVaCL9DGtaAKfv76C0xfao3rft90TTGnfz7wLtUxv20kK3QCY4tx3aRMmMzs/1+8yLlGcmzl8oxEw5/z5H66mpsbV1tb6sm8RkURlZhsG+xOfmkMXEUkSCnQRkSShQBcRSRIKdBGRJKFAFxFJEgp0EZEkoUAXEUkSCnQRkSTh24VFZtYEHBjh24uBEzEsJyjUr8SSjP1Kxj5BcvVrunMuNNAK3wJ9NMysdrArpRKZ+pVYkrFfydgnSN5+9acpFxGRJKFAFxFJEoka6I/5XcAYUb8SSzL2Kxn7BMnbr0sk5By6iIhcLlFH6CIi0o8CXUQkSfgS6GZ2p5ntMrM6M3tkgPVmZt8Lr99sZtcO914zKzKzNWa2J/xc2Gfd34Xb7zKzO5KhX2ZWaWatZvZB+LEqwfr1GTPbZmbdZlbTb3uJ/HkN2K94fV5j1KfvmNnOcPv/NLOCPusS+bMasF/x/G7FnHMurg8gFdgLzAQygE1AVb82K4AXAQNuAN4d7r3APwKPhH9+BPh2+OeqcLtMYEb4/alJ0K9KYGsCf17zgLnAa0BNn20l+uc1WL/G/PMawz7dDqSFf/52En23ButXXL5bY/HwY4S+GKhzztU759qBp4GV/dqsBH7mPOuAAjObMsx7VwI/Df/8U+ATfZY/7Zy76JzbB9SFt5Po/YqXMemXc26Hc27XAPtL6M9riH7Fw1j16RXnXGf4/euA8j7bSuTParB+JSw/Ar0MONTndUN4WSRthnpvqXPuCED4uSSK/cVCvPsFMMPMNprZWjNbNvouDGis+jWa/cVCvPsFY/95xaNPf4k3Eo50f7EQ735BfL5bMZfmwz5tgGX9z50crE0k7x3J/mIh3v06AlQ4506a2SLgN2ZW7ZxrGb7UqOjzurxNUD+vMe2TmX0N6AR+EcX+YiHe/YrXdyvm/Aj0BmBan9flQGOEbTKGeO8xM5vinDsS/lXreBT7i4W49ss5dxG4GP55g5ntBeYAtbHpzrA1R9JmqH6NZn+xENd+xenzGrM+mdmfAx8DPuKc6wnEhP+sBupXHL9bsRfvSXu8f0Tq8Q6i9BykqO7X5o+59ADHe8O9F/gOlx48/Mfwz9VceuCmnrE5cBPvfoV6+oF3wOcwUJQo/erz3te49OBhQn9eQ/RrzD+vMfx/8E5gOxDqt62E/qyG6Fdcvltj8fBnp94R6d14R5+/Fl72EPBQ+GcDHg2v39Lvi3HZe8PLJwG/B/aEn4v6rPtauP0u4K5k6BfwKWBb+H/Q94G7E6xfn8QbVV0EjgEvJ8nnNWC/4vV5jVGf6vDmoT8IP1YlyWc1YL/i+d2K9UOX/ouIJAldKSoikiQU6CIiSUKBLiKSJBToIiJJQoEuIpIkFOgiIklCgS4ikiT+PwGst4Yr59xJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(xt_prf, sw_prf)\n",
    "plt.plot(x.value.squeeze()/(steps*dt), sw.value)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
